# defines class objects for extEOStable and functions to read in ANEOS input and output
# READ THIS CODE IF YOU WANT TO UNDERSTAND THE OBJECT STRUCTURES

%run import-modules

Platform:  Darwin Kernel Version 24.0.0: Tue Sep 24 23:37:25 PDT 2024; root:xnu-11215.1.12~1/RELEASE_ARM64_T6030
python version:  3.11.9
matplotlib version:  3.9.1
hvplot version:  0.10.0
numpy version:  2.0.1
pandas version:  2.2.2

# THIS CELL DEFINES FUNCTIONS AND CUSTOMIZATIONS
#
plt.rcParams["figure.figsize"] = (20,10)
plt.rcParams["patch.force_edgecolor"] = True
#
# define index function (must be a more efficient way to do this...)
get_indexes = lambda x, xs: [i for (y, i) in zip(xs, range(len(xs))) if x == y]
#
# this prevents some error messages for inline plotting in a jupyter notebook
%matplotlib inline 
#
def make_error_boxes(ax, xdata, ydata, xerror, yerror, facecolor='r',
                     edgecolor='None', alpha=0.5):

    # Loop over data points; create box from errors at each point
    errorboxes = [Rectangle((x - xe[0], y - ye[0]), xe.sum(), ye.sum())
                  for x, y, xe, ye in zip(xdata, ydata, xerror.T, yerror.T)]

    # Create patch collection with specified colour/alpha
    pc = PatchCollection(errorboxes, facecolor=facecolor, alpha=alpha,
                         edgecolor=edgecolor)

    # Add collection to axes
    ax.add_collection(pc)

    # Plot errorbars
    artists = ax.errorbar(xdata, ydata, xerr=xerror, yerr=yerror,
                          fmt='None', ecolor='k')

    return artists

# Remember the versions used to run python code
# copy this cell into your own notebooks
# add additional modules to the list in the f.write statements

# Record the verions information for these calculation for posterity
from platform import python_version
from matplotlib import __version__ as mplversion
from numpy import __version__ as npversion
from scipy import __version__ as scipyversion
from pandas import __version__ as pdversion
from eos_table import __version__ as etabversion

# Get date and time values in ISO 8601 format
# importing datetime and time module
# Calling the isoformat() function over the current time including time zone:
import datetime
import zoneinfo
DateTime_in_ISOFormat = datetime.datetime.now(zoneinfo.ZoneInfo('America/Los_Angeles')).isoformat('T', 'seconds') #

versionfilename = './python-versions.txt'
with open(versionfilename, 'a') as f:
    f.write('---------------------------\n')
    f.write(DateTime_in_ISOFormat+'\n')
    f.write('python '+python_version()+'\n')
    f.write('numpy '+npversion+'\n')
    f.write('scipy '+scipyversion+'\n')
    f.write('matplotlib '+mplversion+'\n')
    f.write('pandas '+pdversion+'\n')
    f.write('eos_table '+etabversion+'\n')
    f.write('---------------------------\n')

# Printing Today's date and time in ISO format
print(f"Version history file {versionfilename}")
print(f"Date and time {DateTime_in_ISOFormat}")
print('python '+python_version())
print('numpy '+npversion)
print('scipy '+scipyversion)
print('matplotlib '+mplversion)
print('pandas '+pdversion)
print('eos_table '+etabversion)

Version history file ./python-versions.txt
Date and time 2024-10-18T08:51:18-07:00
python 3.11.9
numpy 2.0.1
scipy 1.14.0
matplotlib 3.9.1
pandas 2.2.2
eos_table 1.1.4

# -------------------------------------------------------------
# USER INPUTS
#-------------------------------------------------------------
# NOTEBOOK MODE
#-------------------------------------------------------------
#
DEVMODE = 0 # 1=development mode; 0=release mode
#
#-------------------------------------------------------------
# OPTIONS TO REMOVE TENSION REGION AND CONSTRUCT A GADGET2 TABLE
#-------------------------------------------------------------
#
REMOVETENSIONFLAG=1  # 0- keep tension region; 1- remove tension region
MAKEGADGETTABLE=1  # 1=make table; 0=do not make gadget2 table
#
#-------------------------------------------------------------
# USER MUST SET THESE VALUES TO DEFINE THE TABULATED EOS GRIDS
#-------------------------------------------------------------
#
# ====>>>>>> YOU NEED TO MAKE SURE THESE VALUES MATCH ANEOS.INPUT  <<<<=====
MATERIALNAME = 'silica'
MODELNAME = 'Fused-silica-ANEOS-SLVTv0.2G1'
# Header information must all be compatible with float format
MATID = 1.0        # MATID number
DATE = 220227.     # Date as a single 6-digit number YYMMDD
VERSION = 0.1      # ANEOS Parameters Version number
FMN = 30.          # Formula weight in atomic numbers for silica
FMW = 60.083       # Formula molecular weight (g/mol)
# The following define the default initial state for material in the 201 table
R0REF   = 2.20     # g/cm3 *** R0REF is inserted into the density array
K0REF   = 4.0E11   # dynes/cm2
T0REF   = 298.     # K -- *** T0REF is inserted into the temperature array
P0REF   = 1.E6     # dynes/cm2 -- this defines the principal Hugoniot calculated below
#
# generate output files for other ANEOS models
modelflag = 0   # 0 - don't change anything above; 1 - Melosh MAPS 2007
#
if modelflag == 1:
    MODELNAME = 'Silica_Melosh_2007'
    # Header information must all be compatible with float format
    MATID = 1.0        # MATID number
    DATE = 070513.     # Date as a single 6-digit number YYMMDD
    VERSION = 1.0      # ANEOS Parameters Version number
    FMN = 30.          # Formula weight in atomic numbers for silica
    FMW = 60.083       # Formula molecular weight (g/mol)
    # The following define the default initial state for material in the 201 table
    R0REF   = 2.65     # g/cm3 *** R0REF is inserted into the density array
    K0REF   = 3.58E11  # dynes/cm2
    T0REF   = 298.     # K -- *** T0REF is inserted into the temperature array
    P0REF   = 1.E6     # dynes/cm2 -- this defines the principal Hugoniot calculated below
#
#-------------------------------------------------------------
# EXAMPLE EOS GRID ARRAYS
#-------------------------------------------------------------
#
# Life is easier if the grid arrays are even in length (after adding T0REF and R0REF below)
#
# Flag for easily switching between different table grids during model development
# 0: low res for testing
# 1: high res for production
gridflag = 11
#
#------------------------------------------------------
if gridflag == 99: # low res for testing
    gridtvals  = [0, 10000., 20000.,1.e7] # K -- this array defines the edges of the segments along the grid, 0 is changed to 1 below
    gridtstyle = [1,1,0] # 1 for linear, 0 for log in each segment
    gridtnpts  = [30., 30.,20.] # number of grid points in each segment
    # define the density array for the SESAME AND GADGET2 tables
    gridrvals  = [1.e-20, 0.1, 2.4, 3.0, 10., 30.] # g/cm3 -- this array defines the edges of the segments along the grid
    gridrstyle = [0, 1, 1, 1, 1] # 1 for linear, 0 for log in each segment
    gridrnpts  = [20, 10, 6, 30, 20] # number of grid points in each segment
    # define the specific entropy array for the GADGET2 table
    gridsvals  = [1.e-4,.25] # MJ/K/kg -- this array defines the edges of the segments along the grid
    gridsstyle = [0] # 1 for linear, 0 for log in each segment
    gridsnpts  = [80] # number of grid points in each segment
#
if gridflag == 0: # medium resolution grid <1E7 entries in 301 table. This is about 250,000 values.
    # define the temperature array for the SESAME table
    # 0 K is changed to 1 K
    # 298 K is inserted into the table if it is not present in the original array
    gridtvals  = [0, 10000., 20000.,1.e7] # K -- this array defines the edges of the segments along the grid, 0 is changed to 1 below
    gridtstyle = [1,1,0] # 1 for linear, 0 for log in each segment
    gridtnpts  = [200.,200.,151.] # number of grid points in each segment
    # define the density array for the SESAME AND GADGET2 tables
    # 0 g/cm3 is removed if present
    gridrvals  = [1.e-20, 0.1, 2.4, 3.0, 10., 30.] # g/cm3 -- this array defines the edges of the segments along the grid
    gridrstyle = [0, 1, 1, 1, 1] # 1 for linear, 0 for log in each segment
    gridrnpts  = [150, 76, 55, 206, 40] # number of grid points in each segment
    # define the specific entropy array for the GADGET2 table
    # 0 MJ/K/kg is removed if present
    gridsvals  = [0.,0.03] # MJ/K/kg -- this array defines the edges of the segments along the grid
    gridsstyle = [1] # 1 for linear, 0 for log in each segment
    gridsnpts  = [500] # number of grid points in each segment
#
if gridflag == 1: # new form for the high resolution table entries in 301 table. Try for 1e6 values
    # define the temperature array for the SESAME table
    # 0 K is changed to 1 K
    # T0REF is inserted into the table if it is not present in the original array
    # need to have extra temperature resolution around the triple point
    gridtvals  = [0, 100, 1600, 2000, 10000., 20000.,1.e7] # K -- this array defines the edges of the segments along the grid, 0 is changed to 1 below
    gridtstyle = [1,1,1,1,1,0] # 1 for linear, 0 for log in each segment
    gridtnpts  = [2, 151., 140, 300.,100.,50.] # number of grid points in each segment
    # define the density array for the SESAME AND GADGET2 tables
    # R0REF is inserted into the table if it is not present in the original array
    # need to have extra density resolution around the triple point
    # 0 g/cm3 is removed if present
    gridrvals  = [1.e-20, 0.1, 2.4, 3.0, 10., 30.] # g/cm3 -- this array defines the edges of the segments along the grid
    gridrstyle = [0, 1, 1, 1, 1] # 1 for linear, 0 for log in each segment
    gridrnpts  = [200, 100, 160, 300, 100] # number of grid points in each segment
    # define the specific entropy array for the GADGET2 table
    # 0 MJ/K/kg is removed if present
    gridsvals  = [1.e-4,.25] # MJ/K/kg -- this array defines the edges of the segments along the grid
    gridsstyle = [0] # 1 for linear, 0 for log in each segment
    gridsnpts  = [800] # number of grid points in each segment

if gridflag == 11: # new form for the high resolution table entries in 301 table. Try for 1e6 values
    # define the temperature array for the SESAME table
    # 0 K is changed to 1 K
    # T0REF is inserted into the table if it is not present in the original array
    # need to have extra temperature resolution around the triple point
    gridtvals  = [0, 100, 1600, 2000, 20000.,1.e7] # K -- this array defines the edges of the segments along the grid, 0 is changed to 1 below
    gridtstyle = [1,1,1,1,0] # 1 for linear, 0 for log in each segment
    gridtnpts  = [5, 201., 200, 200.,50.] # number of grid points in each segment
    # define the density array for the SESAME AND GADGET2 tables
    # R0REF is inserted into the table if it is not present in the original array
    # need to have extra density resolution around the triple point
    # 0 g/cm3 is removed if present
    gridrvals  = [1.e-30, 0.1, 2.4, 3.0, 10., 30.] # g/cm3 -- this array defines the edges of the segments along the grid
    gridrstyle = [0, 1, 1, 1, 1] # 1 for linear, 0 for log in each segment
    gridrnpts  = [300, 100, 160, 200, 100] # number of grid points in each segment
    # define the specific entropy array for the GADGET2 table
    # 0 MJ/K/kg is removed if present
    gridsvals  = [1.e-4,.25] # MJ/K/kg -- this array defines the edges of the segments along the grid
    gridsstyle = [0] # 1 for linear, 0 for log in each segment
    gridsnpts  = [800] # number of grid points in each segment
#----------------------------------------
# END USER INPUT
#----------------------------------------

# THIS CELL HAS CODE TO GENERATE THE GRIDS AND WRITES TO A FILE TO BE READ BY ANEOS
print('GENERATING ARRAYS FOR THE INDICES OF THE TABULATED EOS.')
#-------------- the following code generates the temperature array
gridtarr = []
for iir in  range(1,len(gridtvals)):
    if gridtstyle[iir-1] == 1:
        # linear region
        temptarr = np.arange(gridtnpts[iir-1])/gridtnpts[iir-1]*(gridtvals[iir]-gridtvals[iir-1])+gridtvals[iir-1]
    if gridtstyle[iir-1] == 0:
        # log region
        logstep = (np.log10(gridtvals[iir])-np.log10(gridtvals[iir-1]))/gridtnpts[iir-1]
        temptarr = np.power(10.,np.log10(gridtvals[iir-1])+(np.arange(gridtnpts[iir-1]))*logstep)
    gridtarr = np.concatenate((gridtarr,temptarr))
# check if T0REF already in the grid
#indexT0REF = get_indexes(T0REF,gridtarr)
indexT0REF = np.where(abs(T0REF-gridtarr)<1.E-4)[0]
if len(indexT0REF) == 0:
    # insert T0REF K into the array
    indexT0REF = np.where((gridtarr-T0REF)>0.)[0]
    gridtarr = np.concatenate((gridtarr[0:indexT0REF[0]],[T0REF],gridtarr[indexT0REF[0]::]))
    if indexT0REF[0]>1:
        print('Inserted the REFERENCE TEMPERATURE INTO THE ARRAY:',gridtarr[indexT0REF[0]-2:indexT0REF[0]+3])
    else:
        print('Inserted the REFERENCE TEMPERATURE INTO THE ARRAY:',gridtarr[indexT0REF[0]-1:indexT0REF[0]+3])
else:
    print('T0REF already in temperature array: ',T0REF,gridrarr[indexT0REF[0]-2:indexT0REF[0]+2])
# Remove 0 K or convert to 1 K
indexzero = get_indexes(0.,gridtarr)
indexone = get_indexes(1.,gridtarr)
if len(indexzero) == 1 and len(indexone) == 0:
    gridtarr[indexzero] = 1.
    print('Zero kelvin changed to 1 kelvin.')
else:
    gridtarr = np.delete(gridtarr, indexzero)
    print('Deleted zero K index.')
#
#gridtsize = len(grdtarr)
print('Number of temperature points: ',len(gridtarr))
#print(gridtarr)
#------------ the following code generates the density array
gridrarr = []
for iir in range(1,len(gridrvals)):
    if gridrstyle[iir-1] == 1:
        # linear region
        temprarr = np.arange(gridrnpts[iir-1])/gridrnpts[iir-1]*(gridrvals[iir]-gridrvals[iir-1])+gridrvals[iir-1]
    if gridrstyle[iir-1] == 0:
        # log region
        logstep = (np.log10(gridrvals[iir])-np.log10(gridrvals[iir-1]))/gridrnpts[iir-1]
        temprarr = np.power(10.,np.log10(gridrvals[iir-1])+(np.arange(gridrnpts[iir-1]))*logstep)
    gridrarr = np.concatenate((gridrarr,temprarr))
# check if R0REF already in the grid
#indexR0REF = get_indexes(R0REF,gridrarr)
indexR0REF = np.where(abs(R0REF-gridrarr)<1.E-4)[0]
if len(indexR0REF) == 0:
    # insert R0REF  into the array
    indexR0REF = np.where((gridrarr-R0REF)>0.)[0]
    gridrarr = np.concatenate((gridrarr[0:indexR0REF[0]],[R0REF],gridrarr[indexR0REF[0]::]))
    print('Inserted the REFERENCE DENSITY INTO THE ARRAY:',gridrarr[indexR0REF[0]-2:indexR0REF[0]+3])
    #print('NOT INSERTING REFERENCE DENSITY INTO THE ARRAY.')
else:
    print('R0REF already in density array: ',R0REF,gridrarr[indexR0REF[0]-2:indexR0REF[0]+2])
# Remove 0 g/cm3 if present
indexzero = get_indexes(0.,gridrarr)
if len(indexzero) == 1:
    gridrarr = np.delete(gridrarr, indexzero)
    print('Deleted zero g/cm3 index.')
#
#gridrsize = len(gridrarr)
print('Number of density points: ',len(gridrarr))
#print(gridrarr)
#------------ the following code generates the specific entropy array for gadget tables
gridsarr = []
for iir in range(1,len(gridsvals)):
    if gridsstyle[iir-1] == 1:
        # linear region
        tempsarr = np.arange(gridsnpts[iir-1])/gridsnpts[iir-1]*(gridsvals[iir]-gridsvals[iir-1])+gridsvals[iir-1]
    if gridsstyle[iir-1] == 0:
        # log region
        logstep = (np.log10(gridsvals[iir])-np.log10(gridsvals[iir-1]))/gridsnpts[iir-1]
        tempsarr = np.power(10.,np.log10(gridsvals[iir-1])+(np.arange(gridsnpts[iir-1]))*logstep)
    gridsarr = np.concatenate((gridsarr,tempsarr))
# Remove 0 MJ/K/kg if present
indexzero = get_indexes(0.,gridsarr)
if len(indexzero) == 1:
    gridsarr[indexzero]=1.e-20
#    gridsarr = np.delete(gridsarr, indexzero)
    print('replaced zero MJ/K/kg index with 1.E-20')
#gridssize = len(gridsarr)
#
print('Number of specific entropy points: ',len(gridsarr))
#print(gridsarr)

print('Number of words in the standard SESAME 301 table (NR, NT, R, T, P, IE, HFE)=',2.+len(gridrarr)+len(gridtarr)+len(gridrarr)*len(gridtarr)*3.)
#-----------------------------------------------------------
#
if DEVMODE:
    print('WRITING tablegrid.txt FOR ANEOS in the local directory.\n')
    # write the table arrays to a file for the ANEOS driver to read in
    tablegridfile = open("tablegrid.txt","w")  
    tablegridfile.write("{:.6e}".format(MATID)+'\n')  
    tablegridfile.write("{:.6e}".format(DATE)+'\n')  
    tablegridfile.write("{:.6e}".format(VERSION)+'\n')  
    tablegridfile.write("{:.6e}".format(FMN)+'\n')  
    tablegridfile.write("{:.6e}".format(FMW)+'\n')  
    tablegridfile.write("{:.6e}".format(R0REF)+'\n')  
    tablegridfile.write("{:.6e}".format(K0REF)+'\n')  
    tablegridfile.write("{:.6e}".format(T0REF)+'\n')  
    tablegridfile.write("{:.6e}".format(len(gridrarr))+'\n')  # number of density points
    tablegridfile.write("{:.6e}".format(len(gridtarr))+'\n')  # number of temperature points
    for i in range(0,len(gridrarr)):
        tablegridfile.write("{:.6e}".format(gridrarr[i])+'\n')
    for i in range(0,len(gridtarr)):
        tablegridfile.write("{:.6e}".format(gridtarr[i])+'\n')
    tablegridfile.close()
#
#-----------------------------------------------------
# Output GADGET2 Table information
if MAKEGADGETTABLE:
    print('Number of words in the extended GADGET2 EOS table (NR, NS, R, S, P, T, IE, CS, HFE, CV, KPA, MDQ)=',2.+len(gridrarr)+len(gridtarr)+len(gridrarr)*len(gridtarr)*8.)
#------------------------------------------------------

#----------------------------------------------------
# Sp. Internal Energy grid for Tillotson
#tillearr = np.power(10.,np.arange(600)/600.*6.-3.) # MJ/kg
#tillearr = np.append(0.,tillearr)
#print('Number of specific internal energy points for Tillotson EOS: ',len(tillearr))
#----------------------------------------------------
# Paranoia check
# make sure all table grid values are unique
if (len(np.unique(gridtarr))!=len(gridtarr)):
    print('NON-UNIQUE VALUES IN TEMP ARRAY')
    stopthecode
if (len(np.unique(gridrarr))!=len(gridrarr)):
    print('NON-UNIQUE VALUES IN DENSITY ARRAY')
    stopthecode
if (len(np.unique(gridsarr))!=len(gridsarr)):
    print('NON-UNIQUE VALUES IN SP. ENTROPY ARRAY')
    stopthecode
#if (len(np.unique(tillearr))!=len(tillearr)):
#    print('NON-UNIQUE VALUES IN SP. ENERGY ARRAY')
#    stopthecode

#----------------------------------------------------

print('\nDone generating EOS grids.\n')
print('\nCHECK THAT THE GRID ARRAYS ARE EVEN IN LENGTH.\n')

# END TABLE GRIDDING INFORMATION

GENERATING ARRAYS FOR THE INDICES OF THE TABULATED EOS.
Inserted the REFERENCE TEMPERATURE INTO THE ARRAY: [286.56716418 294.02985075 298.         301.49253731 308.95522388]
Zero kelvin changed to 1 kelvin.
Number of temperature points:  657
Inserted the REFERENCE DENSITY INTO THE ARRAY: [2.17  2.193 2.2   2.216 2.239]
Number of density points:  861
Number of specific entropy points:  800
Number of words in the standard SESAME 301 table (NR, NT, R, T, P, IE, HFE)= 1698551.0
WRITING tablegrid.txt FOR ANEOS in the local directory.

Number of words in the extended GADGET2 EOS table (NR, NS, R, S, P, T, IE, CS, HFE, CV, KPA, MDQ)= 4526936.0

Done generating EOS grids.


CHECK THAT THE GRID ARRAYS ARE EVEN IN LENGTH.

if DEVMODE:
    # call aneos in the local directory
    # output file is in the local directory
    print('RUNNING IN DEVELOPMENT MODE')
    print('Calling aneos in the local directory.')
    completed = subprocess.run(['./aneos-mbam3'])
    print('Aneos call complete. CHECK THAT A NEW ANEOS.OUTPUT FILE WAS CREATED.')
    print('aneos call result: ', completed.returncode)
    print('Check ANEOS.OUTPUT file for errors.')
else:
    print('RUNNING IN RELEASE MODE: NO ANEOS SYSTEM CALL.')
    print('READING IN EOS DISTRIBUTION FILES.')

RUNNING IN DEVELOPMENT MODE
Calling aneos in the local directory.
Aneos call complete. CHECK THAT A NEW ANEOS.OUTPUT FILE WAS CREATED.
aneos call result:  0
Check ANEOS.OUTPUT file for errors.

# READ IN NEW ANEOS MODEL and fill the extEOStable class object
# source in eostable.py
#------------------------------------------------------------------
NewEOS  = etab.extEOStable() # FIRST make new empty EOS object
NewEOS.loadextsesame('NEW-SESAME-EXT.TXT') # LOAD THE EXTENDED 301 SESAME FILE GENERATED BY STSM VERSION OF ANEOS
NewEOS.loadstdsesame('NEW-SESAME-STD.TXT') # LOAD THE STANDARD 301 SESAME FILE GENERATED BY STSM VERSION OF ANEOS
NewEOS.MODELNAME = MODELNAME # string set above in user input
NewEOS.MDQ = np.zeros((NewEOS.NT,NewEOS.ND)) # makes the empty MDQ array
#print(NewEOS.units) # these are the default units for SESAME rho-T tables
#'Units: g/cm3, K, GPa, MJ/kg, MJ/kg, MJ/K/kg, cm/s, MJ/K/kg, KPA flag. 2D arrays are (NT,ND).'

# Add the header info to the table. This could be done during the loading. 
# if made from this notebook, these values are set in the user-input above.
# ** MAKE SURE THEY MATCH ANEOS.INPUT **
NewEOS.MATID   = MATID
NewEOS.DATE    = DATE
NewEOS.VERSION = VERSION
NewEOS.FMN     = FMN
NewEOS.FMW     = FMW
NewEOS.R0REF   = R0REF
NewEOS.K0REF   = K0REF
NewEOS.T0REF   = T0REF
NewEOS.P0REF   = P0REF
#
# Load the information from ANEOS.INPUT and ANEOS.OUTPUT
NewEOS.loadaneos(aneosinfname='ANEOS.INPUT',aneosoutfname='ANEOS.OUTPUT')
#
# calculate the principal Hugoniot from the table
if DEVMODE:
    # calculate the Hugoniot and write to a file. This is from the ANEOS REFERENCE STATE
    NewEOS.calchugoniot(r0=NewEOS.R0REF,t0=NewEOS.T0REF,writefilename='NEW-SESAME-HUG.TXT')
    # calculate the Hugoniot and write to a file pyrolitic Glass initial state, which is not on this EOS surface (proxy crystalline solid + liquid + vapor)
    # this result is basically the same as the requested Hugoniot at r0 and t0, where ANEOS interprets the initial density to be in tension. Differences at very low pressures.
    #NewEOS.calcOffEOSHugoniot(r0=2.20,t0=298.,p0=1.e-4,e0=0.1565,r1=2.65,writefilename='NEW-SESAME-HUG-GLASS.TXT')
    NewEOS.calcOffEOSHugoniot(r0=2.20,t0=298.,p0=1.e-4,e0=0.,r1=2.65,writefilename='NEW-SESAME-HUG-GLASS.TXT')
else:
    # just calculate the Hugoniot and put it in the EOS object
    NewEOS.calchugoniot(r0=NewEOS.R0REF,t0=NewEOS.T0REF)
#    NewEOS.calcOffEOSHugoniot(r0=2.20,t0=298.,p0=1.e-4,e0=0.1565,r1=2.65,writefilename='NEW-SESAME-HUG-GLASS.TXT')
    NewEOS.calcOffEOSHugoniot(r0=2.20,t0=298.,p0=1.e-4,e0=0.,r1=2.65,writefilename='NEW-SESAME-HUG-GLASS.TXT')
#
# 9390/.06/1.e6=0.1565 Richet 9390 J/mol -> MJ/kg
# calculate the 1-bar profile; loop over temp
NewEOS.onebar.T = np.zeros(NewEOS.NT)
NewEOS.onebar.S = np.zeros(NewEOS.NT)
NewEOS.onebar.rho = np.zeros(NewEOS.NT)
it0 = np.where(NewEOS.T >= NewEOS.T0REF)[0]
id0 = np.arange(NewEOS.ND)#np.where(NewEOS.rho >= 0.8*NewEOS.R0REF)[0]
for iit in range(0,NewEOS.NT):
    NewEOS.onebar.T[iit] = NewEOS.T[iit]
    NewEOS.onebar.S[iit] = np.interp(1.E-4,NewEOS.P[iit,id0],NewEOS.S[iit,id0])
    NewEOS.onebar.rho[iit] = np.interp(1.E-4,NewEOS.P[iit,id0],NewEOS.rho[id0])
    #print(iit,NewEOS.onebar.rho[iit],NewEOS.onebar.T[iit],NewEOS.onebar.S[iit]*1.E3)

Done loading ANEOS files.
ANEOS WAS CALLED WITH THE FOLLOWING INPUT, LOADED FROM FILE  ANEOS.INPUT
 *         EOS#      string for EOS name           .....Hugoniot initial........ long
 *                                                 T(eV)           Rho(g/cc)     flag
 * Reference solid Hugoniot at STP initial density and temperature
 ANEOS1    -1        'FUSED_SILICA_v0.2'           THUG=-1         RHUG=2.2      LONG

 *
 *         nelem     type      rho0      temp0(K)  press0    B0/-Cbulk Gamma0    Tdebye(K)
 *         V01       V02       V03       V04       V05       V06       V07       V08
 ANEOS2    2         4         2.20      298.      1.e6      0.40E12   0.57      -1600.

rho0= 2.2
 *         TG model  Gamma     Esep      Tmelt(K)  ..Critical Pt..     Thermal Conduction
 *         or S1     3*C24     erg/g     or -Em    C53       C54       H0        C41
 *         V09       V10       V11       V12       V13       V14       V15       V16
 ANEOS3    1.        6.0       1.20E11   1996.     -10.0E11  0.90      0.        0.

 *         solid     .....Solid-solid phase transition data.....                 rhol/rhos
 *         rhomin    D1        D2        D3        D4        D5        Hfusion   or -rholiq
 *         V17       V18       V19       V20       V21       V22       V23       V24
 ANEOS4    0.        2.92      4.30      0.4E11    2.85E12   1.E1      0.66E10   0.93

 *         Cold comp. limits    ...Liquid EOS params...       ...Interpolation params...
 *         Up        L0        alpha     beta      gamma     C60       C61       C62
 *         V25       V26       V27       V28       V29       V30       V31       V32
 ANEOS5    0.        0.        0.        0.        0.        0.0       -0.851     0.5

 *         Ioniz.    Reactive_chemistry  ...............Molecular Clusters...............
 *         model     Eshift    Sshift    Atoms     Ebind(eV) RotDOF    Rbond     VibDOF
 *         V33       V34       V35       V36       V37       V38       V39       V40
 ANEOS6    0         0.        0.        2         5.0       2.0       1.5E-8    1

 *         Theta_vib Pcold/Expanded      fcv*3nR   QCC1      QCC6      Unused    Unused
 *         (K)       Mieflag   a_exp     adjust    rhomin    psitrans  
 *         V41       V42       V43       V44       V45       V46       V47       V48
 ANEOS7    2000.     1         1.70      1.29      1.E-30    1.D5      0.        0.

 *
 *         Elements in material (number must = nelem)
 *         Z         Atomic Fraction
 ANEOS8    8         0.66667

 ANEOS10   14        0.33333

READING DATA FROM ANEOS OUTPUT FILE  ANEOS.OUTPUT
Writing Hugoniot to file:  NEW-SESAME-HUG.TXT
r1= 2.65
468 32 2.2 298.0 0.0001 0.0 0.0 1.4440455300000001
Writing Hugoniot to file:  NEW-SESAME-HUG-GLASS.TXT

# Make a NOTENSION TABLE and GADGET TABLE
################### REMOVE TENSION REGION ###########################
## The GADGET2 tables do not have a tension region. 
## If the gadget2 table is selected, then must generate a no-tension table
if MAKEGADGETTABLE:
    REMOVETENSIONFLAG=1  # 0- keep tension region; 1- remove tension region

# minimum pressure 1.E-30
# minimum sp energy about 1.E-4 before negative entropy region
# minimum sp entropy about 1.E-5 before negative entropy region
# minimum HFE about 1.E-3 before negative entropy region
# minimum sound speed about 1.E-10 before negative entropy region
# minimum sp heat capacity about 1.E-4 before negative entropy region

## FUTURE WORK: THIS SECTION COULD USE THE LEVER RULE TO FILL IN THE MISSING POINTS
## SINCE MANY VALUES ARE NEAR THE LOWER LIMITS IN THE SUBLIMATION REGION ANYWAY, 
## HERE, TENSION REGION VALUES ARE FILLED IN WITH THE ADJACENT S-V DOME VALUES
## This is obviously bad for interpolations. Fix for future work.
if REMOVETENSIONFLAG:
    # First make a deep copy of the original table with the tension region
    NewEOSorig = deepcopy(NewEOS)
    NewEOS.MODELNAME = NewEOS.MODELNAME+'notension'
    # THIS SECTION REWRITES NEGATIVE PRESSURE REGION
    # The tension region is in the solid-vapor phase boundary. 
    # Use the liquid-vapor phase boundary. Below the temperature of the triple point, it is the sublimation curve.
    # But the pressures for the solid calculation is spurious at very low pressures, so use the vapor pressure.
    # HFE = GFE + TS
    print('REMOVING TENSION REGION FROM THE ORIGINAL ANEOS TABLE.')
    for itt in range(NewEOS.NT-1,-1,-1):
        tmp = np.where(NewEOS.P[itt,:] < 0.)
        if (len(tmp[0]) > 0):
            badind = tmp[0]
            if badind[0] > 0:
                #print(itt, len(badind), NewEOS.T[itt], 'tension only', NewEOS.P[itt,badind[0]-4:badind[0]+1])
                # at this temperature, all densities are the same pressure in the vapor dome
                # so use the adjacent positive pressure value
                # do the same for all the variables since there is very little change in this region of the dome.
                # Errors introduced by interpolating the vapor curve are comparable to
                # recalculating the lever rule. Would need to call ANEOS again to get good end point values
                # to calculate the mass-weighted value via the lever rule.
                NewEOS.P[itt,badind] = np.full(len(badind),NewEOS.P[itt,badind[0]-1])
                NewEOS.U[itt,badind] = np.full(len(badind),NewEOS.U[itt,badind[0]-1])
                NewEOS.S[itt,badind] = np.full(len(badind),NewEOS.S[itt,badind[0]-1])
                NewEOS.A[itt,badind] = np.full(len(badind),NewEOS.A[itt,badind[0]-1])
                NewEOS.cs[itt,badind] = np.full(len(badind),NewEOS.cs[itt,badind[0]-1])
                NewEOS.cv[itt,badind] = np.full(len(badind),NewEOS.cv[itt,badind[0]-1])
                NewEOS.KPA[itt,badind] = np.full(len(badind),2)
                NewEOS.MDQ[itt,badind] = np.full(len(badind),-1) # replaced tension region flag
            else:
                print('SOMETHING VERY WRONG WITH TENSION IN THE TABLE')
                makethecodestop # please tell me a more graceful way to stop the notebook
    #
    # THIS SECTION RESETS NEGATIVE SPECIFIC ENTROPIES WITH ADJACENT VALUES OR MINIMUM POSITIVE VALUES
    # THESE VALUES TYPICALLY DEFINE THE COLD COMPRESSION LIMIT OF THE MATERIAL
    # OR SPURIOUS NEGATIVE ENTROPIES AT VERY LOW PRESSURES (IMPERFECT DEBYE MODEL IN ANEOS)
    for itt in range(NewEOS.NT-1,-1,-1):
        tmp = np.where(NewEOS.S[itt,:] < 0.)
        if (len(tmp[0]) > 0):
            badind = tmp[0]
            if badind[0] > 0:
                #print(itt, antarr[itt], andarr[badind[0]], badind[0],len(badind),'negative entropy', ansarr[itt,badind[0]],ansarr[itt,badind[0]-1])
                NewEOS.P[itt,badind] = np.full(len(badind),NewEOS.P[itt,badind[0]-1]) # use the adjacent pressure at this temperature in the dome
                NewEOS.U[itt,badind] = np.full(len(badind),NewEOS.U[itt,badind[0]-1])
                NewEOS.S[itt,badind] = np.full(len(badind),NewEOS.S[itt,badind[0]-1])
                NewEOS.A[itt,badind] = np.full(len(badind),NewEOS.A[itt,badind[0]-1])
                NewEOS.cs[itt,badind] = np.full(len(badind),NewEOS.cs[itt,badind[0]-1])
                NewEOS.cv[itt,badind] = np.full(len(badind),NewEOS.cv[itt,badind[0]-1])
                #NewEOS.KPA[itt,badind]  # leave as original region flag
                NewEOS.MDQ[itt,badind] = np.full(len(badind),-2) # negative entropies flag
                #print('Some NEGATIVE  S AT THIS INDEX, TEMP(K): ',itt, NewEOS.T[itt],NewEOS.P[itt,badind[0]-1],NewEOS.P[itt,badind[0]-2]) #, badind[0],len(badind),'negative entropy', NewEOS.S[itt,badind[0]])
            if badind[0] == 0:
                print('All S VALUES NEGATIVE AT THIS INDEX, TEMP(K): ',itt, NewEOS.T[itt]) #, badind[0],len(badind),'negative entropy', NewEOS.S[itt,badind[0]])
                NewEOS.S[itt,badind] = np.full(len(badind),1.E-10)            
                NewEOS.U[itt,badind] = np.full(len(badind),1.E-10) # also replace the negative sp. internal energies           
                NewEOS.MDQ[itt,badind] = np.full(len(badind),-2) # negative entropies flag

    #### ONLY WRITE A NEW TABLE IN DEVELOPMENT MODE
    if DEVMODE:
        #### WRITE NEW SESAME TABLE WITH NO TENSION REGION
        #     WRITE STANDARD (SHORT) SESAME FILE
        #     WRITE SESAME TABLE TO FILE LIMITED TO P, E, HFE
        NewEOS.writestdsesame(writestdsesfname="NEW-SESAME-STD-NOTENSION.TXT")
        NewEOS.writeextsesame(writeextsesfname="NEW-SESAME-EXT-NOTENSION.TXT")
    if DEVMODE:
        #### WRITE NEW SESAME TABLE WITH NO TENSION REGION
        #     WRITE STANDARD (SHORT) SESAME FILE
        #     WRITE SESAME TABLE TO FILE LIMITED TO P, E, S
        NewEOS.writegadgetinitsesame(writegadgetinitsesfname="NEW-SESAME-GADGETINIT-NOTENSION.TXT")

#### ONLY WRITE A NEW TABLE IN DEVELOPMENT MODE
if DEVMODE:
    #### WRITE MDQ Flag as a 301-style table
    NewEOS.writemdqsesame(writemdqsesfname="NEW-SESAME-MDQ.TXT")
#
#--------------------------------------------------------------------------
########## MAKE GADGET2 EOS TABLE ARRAYS ###################
# LATER WE CAN TURN THIS INTO A FUNCTION IN THE GADGET EOS OBJECT
# FLAG IS SET ABOVE IN THE USER INPUT SECTION
# MAKEGADGETTABLE = 1 # 1 - do it ; 0 - skip it
if MAKEGADGETTABLE:
    print('Starting interpolation for GADGET table.')
    GadEOS  = etab.extGADtable() # make new empty Gadget EOS object
    #
    # rho array is andarr to match the SESAME table
    # S array is gridsarr
    GadEOS.ND  = len(gridrarr)
    GadEOS.NS  = len(gridsarr)
    GadEOS.rho = gridrarr
    GadEOS.S   = gridsarr
    GadEOS.P   = np.zeros((GadEOS.NS,GadEOS.ND))
    GadEOS.T   = np.zeros((GadEOS.NS,GadEOS.ND))
    GadEOS.U   = np.zeros((GadEOS.NS,GadEOS.ND))
    GadEOS.A   = np.zeros((GadEOS.NS,GadEOS.ND))
    GadEOS.cs  = np.zeros((GadEOS.NS,GadEOS.ND))
    GadEOS.cv  = np.zeros((GadEOS.NS,GadEOS.ND))
    GadEOS.KPA = np.zeros((GadEOS.NS,GadEOS.ND))
    GadEOS.MDQ = np.zeros((GadEOS.NS,GadEOS.ND))
    GadEOS.MODELNAME = NewEOS.MODELNAME+'-gadget-linear-int'
    #
    # some combinations of density and sp. entropy are not in the EOS surface
    # in these cases, apply minimum values
    gadpmin  = np.min(NewEOS.P) # GPa
    gadsmin  = 1.E-30     # MJ/K/kg
    gademin  = np.min(NewEOS.U)    # MJ/kg
    gadtmin  = 1.     # K
    gadcmin  = np.min(NewEOS.cs) # cm/s
    gadcvmin = 0.     # MJ/K/kg
    gadamin  = np.min(NewEOS.A) # can be negative    # MJ/kg
    print('Minimum values used in the rational function interpolation scheme:')
    print('gadpmin,gadsmin,gademin,gadtmin,gadcmin,gadcvmin,gadamin=')
    print(gadpmin,gadsmin,gademin,gadtmin,gadcmin,gadcvmin,gadamin)
    #
    # some combinations of density and sp. entropy are not in the EOS surface
    # in these cases, fill the region with the adjacent real values from the original T-RHO table
    # ----------------------------------------------------------
    # SIMPLE LINEAR INTERPOLATION TO FLIP T-RHO to S-RHO table
    print("LINEAR INTERPOLATION FOR GADGET TABLE.")
    for iid in range(0,GadEOS.ND):
        # same density arrays for the T-rho table and S-rho table
        sall1d = NewEOS.S[:,iid]
        pall1d = NewEOS.P[:,iid]
        eall1d = NewEOS.U[:,iid]
        tall1d = NewEOS.T[:]
        call1d = NewEOS.cs[:,iid]
        cvall1d = NewEOS.cv[:,iid]
        aall1d = NewEOS.A[:,iid]
        kall1d = NewEOS.KPA[:,iid]
        mdqall1d = NewEOS.MDQ[:,iid]
        # interpolate P at the new S
        testinterp = interpolate.interp1d(sall1d,pall1d,bounds_error=False)
        tmpind = np.where(gridsarr < np.max(sall1d))[0]
        GadEOS.P[tmpind,iid] = testinterp(gridsarr[tmpind])
        tmpind = np.where(gridsarr >= np.max(sall1d))[0]
        GadEOS.P[tmpind,iid] = pall1d[NewEOS.NT-1]
        # interpolate E at the new S
        testinterp = interpolate.interp1d(sall1d,eall1d,bounds_error=False)
        tmpind = np.where(gridsarr < np.max(sall1d))[0]
        GadEOS.U[tmpind,iid] = testinterp(gridsarr[tmpind])
        tmpind = np.where(gridsarr >= np.max(sall1d))[0]
        GadEOS.U[tmpind,iid] = eall1d[NewEOS.NT-1]
        # interpolate T at the new S
        testinterp = interpolate.interp1d(sall1d,tall1d,bounds_error=False)
        tmpind = np.where(gridsarr < np.max(sall1d))[0]
        GadEOS.T[tmpind,iid] = testinterp(gridsarr[tmpind])
        tmpind = np.where(gridsarr >= np.max(sall1d))[0]
        GadEOS.T[tmpind,iid] = tall1d[NewEOS.NT-1]
        # interpolate sound speed at the new S
        testinterp = interpolate.interp1d(sall1d,call1d,bounds_error=False)
        tmpind = np.where(gridsarr < np.max(sall1d))[0]
        GadEOS.cs[tmpind,iid] = testinterp(gridsarr[tmpind])
        tmpind = np.where(gridsarr >= np.max(sall1d))[0]
        GadEOS.cs[tmpind,iid] = call1d[NewEOS.NT-1]
        # interpolate sp. heat capacity at the new S
        testinterp = interpolate.interp1d(sall1d,cvall1d,bounds_error=False)
        tmpind = np.where(gridsarr < np.max(sall1d))[0]
        GadEOS.cv[tmpind,iid] = testinterp(gridsarr[tmpind])
        tmpind = np.where(gridsarr >= np.max(sall1d))[0]
        GadEOS.cv[tmpind,iid] = call1d[NewEOS.NT-1]
        # interpolate HFE at the new S
        testinterp = interpolate.interp1d(sall1d,aall1d,bounds_error=False)
        tmpind = np.where(gridsarr < np.max(sall1d))[0]
        GadEOS.A[tmpind,iid] = testinterp(gridsarr[tmpind])
        tmpind = np.where(gridsarr >= np.max(sall1d))[0]
        GadEOS.A[tmpind,iid] = call1d[NewEOS.NT-1]
        # interpolate phase information KPA flag at the new S
        testinterp = interpolate.interp1d(sall1d,kall1d,bounds_error=False)
        tmpind = np.where(gridsarr < np.max(sall1d))[0]
        GadEOS.KPA[tmpind,iid] = testinterp(gridsarr[tmpind])
        tmpind = np.where(gridsarr >= np.max(sall1d))[0]
        GadEOS.KPA[tmpind,iid] = call1d[NewEOS.NT-1]
        # interpolate model development quality flag at the new S
        testinterp = interpolate.interp1d(sall1d,mdqall1d,bounds_error=False)
        tmpind = np.where(gridsarr < np.max(sall1d))[0]
        GadEOS.MDQ[tmpind,iid] = testinterp(gridsarr[tmpind])
        tmpind = np.where(gridsarr >= np.max(sall1d))[0]
        GadEOS.MDQ[tmpind,iid] = call1d[NewEOS.NT-1]
    #
    GadEOSlin = deepcopy(GadEOS)
    #
    '''
    # ----------------------------------------------------------
    # 1D RATIONAL FUNCTION INTERPOLATION METHOD FROM KERLEY
    print("RATIONAL FUNCTION INTERPOLATION TO GENERATE GADGET2 EOS TABLE.")
    for iid in range(0,GadEOS.ND):
        GadEOS.MODELNAME = NewEOS.MODELNAME+'-gadget-ratfunc-int'
        #
        # same density arrays for the T-rho table and S-rho table
        sall1d  = NewEOS.S[:,iid]
        pall1d  = NewEOS.P[:,iid]
        eall1d  = NewEOS.U[:,iid]
        tall1d  = NewEOS.T[:]
        call1d  = NewEOS.cs[:,iid]
        cvall1d = NewEOS.cv[:,iid]
        aall1d  = NewEOS.A[:,iid]
        #kall1d = NewEOS.KPA[:,iid]
        #mdqall1d = NewEOS.MDQ[:,iid]
        #
        # interpolate P at the new S
        GadEOS.P[:,iid] = rf1d(sall1d,pall1d,gridsarr,gadpmin)
        # interpolate E at the new S
        GadEOS.U[:,iid] = rf1d(sall1d,eall1d,gridsarr,gademin)
        # interpolate T at the new S
        GadEOS.T[:,iid] = rf1d(sall1d,tall1d,gridsarr,gadtmin)
        # interpolate sound speed at the new S
        GadEOS.cs[:,iid] = rf1d(sall1d,call1d,gridsarr,gadcmin)
        # interpolate sp. heat capacity at the new S
        GadEOS.cv[:,iid] = rf1d(sall1d,cvall1d,gridsarr,gadcvmin)
        # interpolate HFE at the new S
        GadEOS.A[:,iid] = rf1d(sall1d,aall1d,gridsarr,gadamin)
        # Do not interpolate KPA or MDQ with rational function interp. Use Linear.
        # do not overwrite the linear interpolation above.
    '''
    #
    print('Done interpolating the gadget2 table.')
    #print('Using RATIONAL FUNCTION INTERPOLATION EXCEPT FOR KPA and MDQ FLAGS (bilinear).')
    if DEVMODE:
        # ONLY WRITE GADGET2 TABLE TO DISK IN DEVELOPMENT MODE
        #### WRITE NEW GADGET2 TABLE WITH NO TENSION REGION
        #### This ascii table format follows the sesame table format of 5 columns of numbers in 16.8e format
        #### the STYLE variable is just a counter to 5 to know when to write a newline
        GadEOS.writestdgadget(writestdgadgetfname="NEW-GADGET-STD-NOTENSION.TXT")
        GadEOS.writeextgadget(writeextgadgetfname="NEW-GADGET-EXT-NOTENSION.TXT")

REMOVING TENSION REGION FROM THE ORIGINAL ANEOS TABLE.
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  10 137.3134
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  9 129.8507
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  8 122.3881
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  7 114.9254
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  6 107.4627
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  5 100.0
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  4 80.0
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  3 60.0
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  2 40.0
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  1 20.0
All S VALUES NEGATIVE AT THIS INDEX, TEMP(K):  0 1.0
Done writing the STD SESAME 301 table to local directory:  NEW-SESAME-STD-NOTENSION.TXT
Done writing the EXT SESAME 301 table to local directory:  NEW-SESAME-EXT-NOTENSION.TXT
Done writing the GADGET INIT SESAME 301 notension table to local directory:  NEW-SESAME-GADGETINIT-NOTENSION.TXT
Done writing the MDQ Flag as a 301-style table to local directory:  NEW-SESAME-MDQ.TXT
Starting interpolation for GADGET table.
Minimum values used in the rational function interpolation scheme:
gadpmin,gadsmin,gademin,gadtmin,gadcmin,gadcvmin,gadamin=
9.99999968e-31 1e-30 1e-10 1.0 1e-10 0.0 -3116194.75
LINEAR INTERPOLATION FOR GADGET TABLE.
Done interpolating the gadget2 table.
Done writing the STD GADGET TABLE FILE with notension to local directory:  NEW-GADGET-STD-NOTENSION.TXT
Done writing the EXTENDED GADGET TABLE FILE with notension to local directory:  NEW-GADGET-EXT-NOTENSION.TXT

%%capture 
## This cell loads the modules, data, and EOS parameters
## capture magic function hides the output; comment out to see diagnostic info during data load

# These scripts and data will be distributed in a separate GitHub repository (Amodeo et al., in prep.)

qtz_liq_hug = 'KD09' # or `KD13` # pick which quartz Hugoniot equation
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/load_paths
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/load_modules
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/load_silica_data
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/define_eos
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/define_mc
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/load_vapor_domes.py
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/load_liquid_model.py
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/load_solid_models.py
%run -i /Users/sts/Library/CloudStorage/Dropbox-Personal/eos/github-group/silica-sm/silica-book/scripts/hugoniot-calcs.py

# Read in additional SiO2 system data that is not read in load_silica_data.py
# silicadatadir = SET IN LOAD_SILICA_DATA.PY

# Use these units:
units_str = 'rho g/cm3, T K, P GPa, U MJ/kg, A MJ/kg, S kJ/K/kg, cs cm/s, cv MJ/K/kg, us km/s, up km/s'

# Bourova, E., Richet, P., 1998. Quartz and Cristobalite: high-temperature cell parameters and volumes of fusion. Geophys. Res. Lett. 25, 2333–2336. https://doi.org/10.1029/98GL01581
# Bourova & Richet 1998 volumes of fusion qtz, cristobalite, liquid
qtz_mc_BR98_file = silicadatadir+'BourovaRichet1998-qtz-melt-curve.txt'
qtz_mc_BR98 = etab.EOSmeltcurve()
qtz_mc_BR98.Pl = np.loadtxt(qtz_mc_BR98_file,delimiter=',',skiprows=2,usecols=[0])*1.E8/1.e9 # kbar -> GPa
qtz_mc_BR98.Ps = np.loadtxt(qtz_mc_BR98_file,delimiter=',',skiprows=2,usecols=[0])*1.E8/1.e9 # kbar -> GPa
qtz_mc_BR98.T = np.loadtxt(qtz_mc_BR98_file,delimiter=',',skiprows=2,usecols=[1]) # K
qtz_mc_BR98.units = units_str
qtz_mc_BR98.label = 'Melt curve Bourova & Richet 1998'
#print(qtz_mc_BR98.Pl)

# Millot, M., Dubrovinskaia, N., Černok, A., Blaha, S., Dubrovinsky, L., Braun, D.G., Celliers, P.M., Collins, G.W., Eggert, J.H., Jeanloz, R., 2015. Shock compression of stishovite and melting of silica at planetary interior conditions. Science (80-. ). 347, 418–420. https://doi.org/10.1126/science.1261507
# Millot et al. Science 2015 Melt curve fitted Simon's equation using stishovite crossing point
#stv_mc_MM15_file   = silicadatadir+'Millot-melt-curve.txt'
#stv_mc_MM15 = EOSmeltcurve()
#stv_mc_MM15.Pl = np.loadtxt(stvmcfile,delimiter=',',skiprows=0,usecols=[0])*1.e9 # Pa
#stv_mc_MM15.Ps = stv_mc.Pl
#stv_mc_MM15.T = np.loadtxt(stvmcfile,delimiter=',',skiprows=0,usecols=[1])*1.e3 # K
stv_mc_MM15 = etab.EOSmeltcurve()
stv_mc_MM15.Pl = np.arange(0,600) # GPa
stv_mc_MM15.Ps = np.arange(0,600) # GPa
stv_mc_MM15.T = 1968.5 + 307.8*np.power(stv_mc_MM15.Pl,0.485) # K with P in GPa
stv_mc_MM15.units = units_str
stv_mc_MM15.label = 'Melt curve Millot et al. 2015'
#print(stv_mc_MM15.T)

# McCoy, C.A., Gregor, M.C., Polsin, D.N., Fratanduono, D.E., Celliers, P.M., Boehly, T.R., Meyerhofer, D.D., 2016. Measurements of the sound velocity of shock-compressed liquid silica to 1100 GPa. J. Appl. Phys. 120. https://doi.org/10.1063/1.4972338
# McCoy et al. JAP 2016 Gammas from sound speeds on the Hugoniot
qtz_hug_McCoy16_file = silicadatadir+'McCoy2016-qtz-hug-gamma.txt'
qtz_hug_McCoy16 = etab.EOShugoniot()
qtz_hug_McCoy16.P = np.loadtxt(qtz_hug_McCoy16_file,delimiter=' ',skiprows=2,usecols=[1]) # GPa
qtz_hug_McCoy16.cs = np.loadtxt(qtz_hug_McCoy16_file,delimiter=' ',skiprows=2,usecols=[5]) # km/s
qtz_hug_McCoy16.gamma = np.loadtxt(qtz_hug_McCoy16_file,delimiter=' ',skiprows=1,usecols=[7]) # [-] 
qtz_hug_McCoy16.rho = np.interp(qtz_hug_McCoy16.P,qtz_hug_KD.P,qtz_hug_KD.rho) # g/cm3
qtz_hug_McCoy16.units = units_str
qtz_hug_McCoy16.label = 'Qtz McCoy et al. 2016'
#print(mccoy_hug.S)
#print(mccoy_hug.rho)

# NIST-JANAF web site SiO2
qtz_1bar_JANAF_file = silicadatadir+'JANAF-qtz-crystal-liquid.txt'
qtz_1bar_JANAF = etab.EOS1barcurve()
qtz_1bar_JANAF.T = np.loadtxt(qtz_1bar_JANAF_file,delimiter='\t',skiprows=2,usecols=[0]) # K
qtz_1bar_JANAF.S = np.loadtxt(qtz_1bar_JANAF_file,delimiter='\t',skiprows=2,usecols=[2])/0.0600843/1.e3 # J/K/mol -> kJ/K/kg
qtz_1bar_JANAF.units = units_str
qtz_1bar_JANAF.label = '1 bar JANAF'
#print(qtz_1bar_JANAF.S)
#janafgasfile = silicadatadir+'janaf-gas.txt'
#janafgas_1bar = EOS1barcurve()
#janafgas_1bar.T = np.loadtxt(janafgasfile,delimiter='\t',skiprows=3,usecols=[0]) # K
#janafgas_1bar.S = np.loadtxt(janafgasfile,delimiter='\t',skiprows=3,usecols=[2])/0.0600843 # J/K/mol -> J/K/kg
#print(janafgas_1bar.S) # this does not seem to be right?? does not agree with Melosh 2007

## read in Connolly data
C16highs = np.loadtxt(silicadatadir+'Connolly-high-cp.csv',delimiter=',',usecols=[0])
C16hight = np.loadtxt(silicadatadir+'Connolly-high-cp.csv',delimiter=',',usecols=[1])
C16bsrho  = np.loadtxt(silicadatadir+'Connolly-Fig4b-black-T-rho.csv',delimiter=',',usecols=[0])/1.e3
C16bstemp = np.loadtxt(silicadatadir+'Connolly-Fig4b-black-T-rho.csv',delimiter=',',usecols=[1])


C16lows = np.loadtxt(silicadatadir+'Connolly-low-cp.csv',delimiter=',',usecols=[0])
C16lowt = np.loadtxt(silicadatadir+'Connolly-low-cp.csv',delimiter=',',usecols=[1])
C16bdrho  = np.loadtxt(silicadatadir+'Connolly-Fig4b-black-dashed-T-rho.csv',delimiter=',',usecols=[0])/1.e3
C16bdtemp = np.loadtxt(silicadatadir+'Connolly-Fig4b-black-dashed-T-rho.csv',delimiter=',',usecols=[1])

C16sils = np.loadtxt(silicadatadir+'Connolly-SiL-VWRK.csv',delimiter=',',usecols=[0])
C16silt = np.loadtxt(silicadatadir+'Connolly-SiL-VWRK.csv',delimiter=',',usecols=[1])
C16rsrho  = np.loadtxt(silicadatadir+'Connolly-Fig4b-red-T-rho.csv',delimiter=',',usecols=[0])/1.e3
C16rstemp = np.loadtxt(silicadatadir+'Connolly-Fig4b-red-T-rho.csv',delimiter=',',usecols=[1])

# Load Karki gammas 2007
v0 = 45.8 # A/SiO2
rho0 = sio2mol/v0*1.6605391 # g/cm3

# STS estimated melt curve line for data review draw as dotted above 155 GPa
dat_file   = silicadatadir+'karki-gamma.txt'
dat         = etab.EOSmeltcurve()
dat.matid    = [174]
dat.material = ['Quartz']
dat.formula  = ['SiO2']
dat.rho_ref = [2.650] # g/cm3
sio2mol = 28.085+ 2*15.999
vol=np.loadtxt(dat_file,delimiter=',',skiprows=1,usecols=[0])*v0
dat.rho      = sio2mol/vol*1.6605391 # amu/a^3 -> g/cm3
dat.rho_units = ['g/cm3']
dat.gamma=np.loadtxt(dat_file,delimiter=',',skiprows=1,usecols=[1])

dat_file   = silicadatadir+'karki-cv.txt'
dat.cv=np.loadtxt(dat_file,delimiter=',',skiprows=1,usecols=[1])*3*8.314/(sio2mol/1000) # J/K/kg

#print(dat.gamma)
#print(dat.rho)
#print(dat.cv)
karkigcvorig = deepcopy(dat)
dat=[]
vol=[]

# PLOT COLOR SHADED VALUES RHO-T TABLE
lw=1
sfont = 15
ds = (10,10)
plt.rcParams['pcolor.shading'] ='nearest'

%matplotlib inline
fig, ([ax0, ax1], [ax2,ax3], [ax4,ax5])= plt.subplots(nrows=3,ncols=2,figsize=(14,24))
plt.subplots_adjust(wspace=0.4)
#plt.tight_layout(pad=0.6)

##================ Pressure
im1=ax0.pcolormesh(NewEOS.rho,NewEOS.T,(np.log10(NewEOS.P)+25)/29.,cmap=local_cmaps.parula,vmin=0,vmax=1,edgecolors='face')
ax0.plot(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax0.plot(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax0.plot(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax0.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax0.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax0.set_ylim(0.,20000)
ax0.set_xlim(0,30.)
ax0.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax0.set_ylabel('Temperature (K)',size=sfont)
ax0.set_title('Colors=log$_{10}$(P (GPa))',size=sfont)
ax0.tick_params(labelsize=sfont)

divider0 = make_axes_locatable(ax0)
cax0 = divider0.append_axes('right', size='5%', pad=0.05)
cbar0 = fig.colorbar(im1, cax=cax0, ticks = [],orientation='vertical')
#cbar0.ax.set_yticklabels(['  ','  ','  '])  # vertical colorbar
cbar0.set_label('Log$_{10}$(P) [GPa]',size=sfont)
#cbar0.ax.set_yticklabels(['-25','-11','+4'])  # vertical colorbar

im1=ax1.pcolormesh(NewEOS.rho,NewEOS.T,(np.log10(NewEOS.P)+25)/29.,cmap=local_cmaps.parula,vmin=0,vmax=1,edgecolors='face')
ax1.loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax1.loglog(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax1.loglog(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax1.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax1.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax1.set_ylim(10.,1.E7)
ax1.set_xlim(1.E-21,100.)
ax1.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax1.set_ylabel('Temperature (K)',size=sfont)
ax1.set_title('Colors=log$_{10}$(P (GPa))',size=sfont)
ax1.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax1)
cax = divider.append_axes('right', size='5%', pad=0.05)
cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
cbar.ax.set_yticklabels(['-25','-11', '+4'])  # vertical colorbar
cbar.set_label('Log$_{10}$(P) [GPa]',size=sfont)


#================ Sp. Int. Energy
im1=ax2.pcolormesh(NewEOS.rho,NewEOS.T,(np.log10(NewEOS.U)+4.)/10.,cmap=local_cmaps.parula,vmin=0,vmax=1,edgecolors='face')
ax2.plot(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax2.plot(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax2.plot(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax2.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax2.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax2.set_ylim(0.,20000)
ax2.set_xlim(0,30.)
ax2.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax2.set_ylabel('Temperature (K)',size=sfont)
ax2.set_title('Colors=log$_{10}$(Sp. Int. Energy (MJ/kg))',size=sfont)
ax2.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax2)
cax = divider.append_axes('right', size='5%', pad=0.05)
# does not want to label ticks in this location
cbar = fig.colorbar(im1, cax=cax, ticks = [], orientation='vertical')
#cbar.ax.set_yticklabels(['','',''])  # vertical colorbar
cbar.set_label('Log$_{10}$(E) [MJ/kg]',size=sfont)
#cbar.ax.set_yticklabels(['-4','+1','+6'])  # vertical colorbar


ax3.pcolormesh(NewEOS.rho,NewEOS.T,(np.log10(NewEOS.U)+4.)/10.,cmap=local_cmaps.parula,vmin=0,vmax=1)
ax3.loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax3.loglog(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax3.loglog(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax3.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax3.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax3.set_ylim(10.,1.E7)
ax3.set_xlim(1.E-21,100.)
ax3.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax3.set_ylabel('Temperature (K)',size=sfont)
ax3.set_title('Colors=log$_{10}$(E (MJ/kg))',size=sfont)
ax3.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax3)
cax = divider.append_axes('right', size='5%', pad=0.05)
#cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
#cbar.ax.set_yticklabels(['-6', '-3.5', '-1'])  # vertical colorbar
cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
cbar.ax.set_yticklabels(['-4', '+1','+6'])  # vertical colorbar
cbar.set_label('Log$_{10}$(E) [MJ/kg]',size=sfont)

#================ Sp. Entropy
im1=ax4.pcolormesh(NewEOS.rho,NewEOS.T,(np.log10(NewEOS.S)+4.)/4.,cmap=local_cmaps.parula,vmin=0,vmax=1,edgecolors='face')
ax4.plot(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax4.plot(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax4.plot(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax4.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax4.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax4.set_ylim(0.,20000)
ax4.set_xlim(0,30.)
ax4.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax4.set_ylabel('Temperature (K)',size=sfont)
ax4.set_title('Colors=log$_{10}$(Sp. Entropy (MJ/K/kg))',size=sfont)
ax4.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax4)
cax = divider.append_axes('right', size='5%', pad=0.05)
cbar = fig.colorbar(im1, cax=cax, ticks = [], orientation='vertical')
#cbar.ax.set_yticklabels(['-4', '-2', '0'])  # vertical colorbar
cbar.set_label('Log$_{10}$(S) [MJ/K/kg]',size=sfont)


ax5.pcolormesh(NewEOS.rho,NewEOS.T,(np.log10(NewEOS.S)+4.)/4.,cmap=local_cmaps.parula,vmin=0,vmax=1)
ax5.loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax5.loglog(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax5.loglog(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax5.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax5.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax5.set_ylim(10.,1.E7)
ax5.set_xlim(1.E-21,100.)
ax5.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax5.set_ylabel('Temperature (K)',size=sfont)
ax5.set_title('Colors=log$_{10}$(Sp. Entropy (MJ/K/kg))',size=sfont)
ax5.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax5)
cax = divider.append_axes('right', size='5%', pad=0.05)
#cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
#cbar.ax.set_yticklabels(['-6', '-3.5', '-1'])  # vertical colorbar
cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
cbar.ax.set_yticklabels(['-4', '-2','0'])  # vertical colorbar
cbar.set_label('Log$_{10}$(S) [MJ/K/kg]',size=sfont)


##============= start new plot
fig2, ([ax0, ax1], [ax2,ax3], [ax4,ax5])= plt.subplots(nrows=3,ncols=2,figsize=(14,24))
plt.subplots_adjust(wspace=0.4)

#================ Sound Speed
im1=ax0.pcolormesh(NewEOS.rho,NewEOS.T,(np.log10(NewEOS.cs)+3.)/10.,cmap=local_cmaps.parula,vmin=0,vmax=1,edgecolors='face')
ax0.plot(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax0.plot(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax0.plot(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax0.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax0.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax0.set_ylim(0.,20000)
ax0.set_xlim(0,30.)
ax0.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax0.set_ylabel('Temperature (K)',size=sfont)
ax0.set_title('Colors=log$_{10}$(Sound Speed (cm/s))',size=sfont)
ax0.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax0)
cax = divider.append_axes('right', size='5%', pad=0.05)
cbar = fig2.colorbar(im1, cax=cax, ticks = [], orientation='vertical')
#cbar.ax.set_yticklabels(['-3', '+2', '+7'])  # vertical colorbar
cbar.set_label('Log$_{10}$(Cs) [cm/s]',size=sfont)


ax1.pcolormesh(NewEOS.rho,NewEOS.T,(np.log10(NewEOS.cs)+3.)/10.,cmap=local_cmaps.parula,vmin=0,vmax=1)
ax1.loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax1.loglog(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax1.loglog(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax1.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax1.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax1.set_ylim(10.,1.E7)
ax1.set_xlim(1.E-21,100.)
ax1.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax1.set_ylabel('Temperature (K)',size=sfont)
ax1.set_title('Colors=log$_{10}$(sound speed (cm/s))',size=sfont)
ax1.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax1)
cax = divider.append_axes('right', size='5%', pad=0.05)
#cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
#cbar.ax.set_yticklabels(['-6', '-3.5', '-1'])  # vertical colorbar
cbar = fig2.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
cbar.ax.set_yticklabels(['-3', '+2', '+7'])  # vertical colorbar
cbar.set_label('Log$_{10}$(Cs) [cm/s]',size=sfont)


#================ HFE
amin = np.min(NewEOS.A)
amin = -1000.
alogmax = 1500. #np.max(anaarrodd)-amin #np.max(np.log10(anaarrodd-amin))
#print('amin, log10(amin)',amin,np.log10(-amin))
#print('alogmax=',alogmax)
im1=ax2.pcolormesh(NewEOS.rho,NewEOS.T,(NewEOS.A-amin)/alogmax,cmap=local_cmaps.parula,vmin=0,vmax=1,edgecolors='face')
ax2.plot(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax2.plot(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax2.plot(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax2.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax2.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax2.set_ylim(0.,20000)
ax2.set_xlim(0,30.)
ax2.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax2.set_ylabel('Temperature (K)',size=sfont)
ax2.set_title('Colors=HFE (MJ/kg)',size=sfont)
ax2.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax2)
cax = divider.append_axes('right', size='5%', pad=0.05)
cbar = fig2.colorbar(im1, cax=cax, ticks = [], orientation='vertical')
#cbar.ax.set_yticklabels(['-1000','-500','0','500'])  # vertical colorbar
cbar.set_label('Log$_{10}$(HFE) [MJ/kg]',size=sfont)


ax3.pcolormesh(NewEOS.rho,NewEOS.T,(NewEOS.A-amin)/alogmax,cmap=local_cmaps.parula,vmin=0,vmax=1)
ax3.loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax3.loglog(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax3.loglog(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax3.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax3.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax3.set_ylim(10.,1.E7)
ax3.set_xlim(1.E-21,100.)
ax3.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax3.set_ylabel('Temperature (K)',size=sfont)
ax3.set_title('Colors=HFE (MJ/kg)',size=sfont)
ax3.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax3)
cax = divider.append_axes('right', size='5%', pad=0.05)
#cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
#cbar.ax.set_yticklabels(['-6', '-3.5', '-1'])  # vertical colorbar
cbar = fig2.colorbar(im1, cax=cax, ticks = [0,.333,.6667,  1], orientation='vertical')
cbar.ax.set_yticklabels(['-1000','-500','0','500'])  # vertical colorbar
cbar.set_label('Log$_{10}$(HFE) [MJ/kg]',size=sfont)



#================ KPA flag
cmap2 = plt.get_cmap('plasma', 6)    # 6 discrete colors

im1=ax4.pcolormesh(NewEOS.rho,NewEOS.T,((NewEOS.KPA)-0.5)/6.,cmap=cmap2,vmin=0,vmax=1,edgecolors='face')
#im1=ax4.pcolormesh(NewEOS.rho,NewEOS.T,((NewEOS.KPA))/6.,cmap=local_cmaps.parula,vmin=0,vmax=1,edgecolors='face')
ax4.plot(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax4.plot(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax4.plot(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax4.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax4.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax4.set_ylim(0.,20000)
ax4.set_xlim(0,30.)
ax4.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax4.set_ylabel('Temperature (K)',size=sfont)
ax4.set_title('Colors=KPA flag (1-6)',size=sfont)
ax4.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax4)
cax = divider.append_axes('right', size='5%', pad=0.05)
cbar = fig2.colorbar(im1, cax=cax, ticks = [], orientation='vertical')
#cbar.ax.set_yticklabels(['0','1','2','3','4','5',  '6'])  # vertical colorbar
cbar.set_label('KPA Flag',size=sfont)


ax5.pcolormesh(NewEOS.rho,NewEOS.T,((NewEOS.KPA)-0.5)/6.,cmap=cmap2,vmin=0,vmax=1)
ax5.loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax5.loglog(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax5.loglog(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax5.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax5.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax5.set_ylim(10.,1.E7)
ax5.set_xlim(1.E-21,100.)
ax5.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax5.set_ylabel('Temperature (K)',size=sfont)
ax5.set_title('Colors=KPA flag (1-6)',size=sfont)
ax5.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax5)
cax = divider.append_axes('right', size='5%', pad=0.05)
#cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
#cbar.ax.set_yticklabels(['-6', '-3.5', '-1'])  # vertical colorbar
#cbar = fig.colorbar(im1, cax=cax, ticks = np.arange(13)/12, orientation='vertical')
cbar = fig2.colorbar(im1, cax=cax, ticks = [], orientation='vertical')
#cbar.ax.set_yticklabels(['','1','','2','','3','','4','','5','',  '6'])  # vertical colorbar
cbar.set_label('KPA Flag',size=sfont)

# check the user-defined model data quality flag across the EOS surface
uniqmdqarr = np.unique(NewEOS.MDQ)
print("Unique values in the MDQ flag array:", uniqmdqarr)
print("Value: # elements with this value")
#get_indexes = lambda x, xs: [i for (y, i) in zip(xs, range(len(xs))) if x == y]

for i in uniqmdqarr:
    ix = get_indexes(i,np.reshape(NewEOS.MDQ,NewEOS.ND*NewEOS.NT))
    print(i,': ',len(ix))


##============= start new plot
fig3, ([ax4, ax5])= plt.subplots(nrows=1,ncols=2,figsize=(14,8))
plt.subplots_adjust(wspace=0.4)

#================ MDQ flag
cmap2 = plt.get_cmap('plasma', 6)    # 6 discrete colors

#im1=ax4.pcolormesh(NewEOS.rho,NewEOS.T,((NewEOS.MDQ)+5)/6.,cmap=local_cmaps.parula,vmin=0,vmax=1,edgecolors='face')
im1=ax4.pcolormesh(NewEOS.rho,NewEOS.T,((NewEOS.MDQ)+5)/6.,cmap=cmap2,vmin=0,vmax=1,edgecolors='face')
ax4.plot(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax4.plot(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax4.plot(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax4.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax4.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax4.set_ylim(0.,20000)
ax4.set_xlim(0,30.)
ax4.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax4.set_ylabel('Temperature (K)',size=sfont)
ax4.set_title('Colors=MDQ flag',size=sfont)
ax4.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax4)
cax = divider.append_axes('right', size='5%', pad=0.05)
cbar = fig3.colorbar(im1, cax=cax, ticks = [], orientation='vertical')
#cbar.ax.set_yticklabels(['','-5','','-4','','-3','','-2','','-1','','0','','1'])  # vertical colorbar

cmap2 = plt.get_cmap('plasma', 6)    # 6 discrete colors
ax5.pcolormesh(NewEOS.rho,NewEOS.T,((NewEOS.MDQ)+5)/6.,cmap=cmap2,vmin=0,vmax=1)
#ax5.pcolormesh(NewEOS.rho,NewEOS.T,((NewEOS.MDQ)+5)/6.,cmap=local_cmaps.parula,vmin=0,vmax=1)
ax5.loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax5.loglog(NewEOS.vc.rv,NewEOS.vc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
ax5.loglog(NewEOS.vc.rl[0],NewEOS.vc.T[0],'ko')
if NewEOS.mc.NT > 0:
    ax5.plot(NewEOS.mc.rl,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)
#    ax5.plot(NewEOS.mc.rs,NewEOS.mc.T,color="black",linewidth=lw,linestyle='--',dashes=ds)

ax5.set_ylim(10.,1.E7)
ax5.set_xlim(1.E-21,100.)
ax5.set_xlabel('Density (g/cm$^3$)',size=sfont)
ax5.set_ylabel('Temperature (K)',size=sfont)
ax5.set_title('Colors=MDQ flag',size=sfont)
ax5.tick_params(labelsize=sfont)

divider = make_axes_locatable(ax5)
cax = divider.append_axes('right', size='5%', pad=0.05)
#cbar = fig.colorbar(im1, cax=cax, ticks = [0, 0.5, 1], orientation='vertical')
#cbar.ax.set_yticklabels(['-6', '-3.5', '-1'])  # vertical colorbar
cbar = fig3.colorbar(im1, cax=cax, ticks = np.arange(13)/12, orientation='vertical')
#cbar.ax.set_yticklabels(['','-5','','-4','','-3','','-2','','-1','','0','','1'])  # vertical colorbar
cbar.ax.set_yticklabels(['','-5','','-4','','-3','','-2','','-1','','0',''])  # vertical colorbar

Unique values in the MDQ flag array: [-2. -1.  0.]
Value: # elements with this value
-2.0 :  29077
-1.0 :  5726
0.0 :  530874

[Text(1, 0.0, ''),
 Text(1, 0.08333333333333333, '-5'),
 Text(1, 0.16666666666666666, ''),
 Text(1, 0.25, '-4'),
 Text(1, 0.3333333333333333, ''),
 Text(1, 0.4166666666666667, '-3'),
 Text(1, 0.5, ''),
 Text(1, 0.5833333333333334, '-2'),
 Text(1, 0.6666666666666666, ''),
 Text(1, 0.75, '-1'),
 Text(1, 0.8333333333333334, ''),
 Text(1, 0.9166666666666666, '0'),
 Text(1, 1.0, '')]

# Display the initial conditions for the Hugoniot in ANEOS.OUTPUT    
print('\nRequested ANEOS Hugoniot initial state from (THUG,RHUG) (red curves in plots below):')
print('Initial Pressure (GPa):',NewEOS.anhug.P[0])
print('Initial Temperature (K):',NewEOS.anhug.T[0])
print('Initial Sp. Entropy (MJ/K/kg):',NewEOS.anhug.S[0])
print('Initial Sp. Internal Energy (MJ/kg):',NewEOS.anhug.U[0])
print('Initial Density (g/cm3):',NewEOS.anhug.rho[0])
#

#========== PLOTTING STYLES FOR LOTS OF PLOTS =================
styleqmd = ['^','red',5,'Fo QMD calc (Root et al. 2018)']
stylezdat = ['^','black',5,'Fo Z data (Root et al. 2018)']
stylegg = ['+','green',5,'Fo Gas gun data\n(Mosenfelder et al. 2007)']
styleliq = ['^','orange',7,'Liquid Fo data corrected\n(Thomas & Asimow 2013)']
styleliqhug = ['-','orange',1,'ANEOS Liquid 2273K Hugoniot']
stylezfit = ['-','green',1,'Fo Z fit (Davies et al. 2019)']
styleanhug = ['-','red',1,'ANEOS Glass-tension Hugoniot']
stylestphug = ['-','blue',1,'ANEOS Qtz Hugoniot']
styleanphase = ['-','black',1,'ANEOS Phase Boundaries']
stylephasept = ['o','black',6]



#======================LOTS OF PLOTS============================
# plot the Hugoniots for comparison
sfont = 15
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(14,18))
plt.subplots_adjust(wspace=0.4)

#-------------
# Pressure-Density Data vs. ANEOS
ai=0
aj=0
#axes[ai,aj].plot(NewEOS.anhug.rho,NewEOS.anhug.P,styleanhug[0],linewidth=styleanhug[2],color=styleanhug[1],label=styleanhug[3])

#axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
#axes[ai,aj].plot(NewEOS.hugo.rho,NewEOS.hugo.P,'-.',color='orange',label='ANEOS FS Hugoniot')

if 0:
    axes[ai,aj].scatter(qtz_hug_W62.rho,qtz_hug_W62.P,marker='+',color='purple',label=qtz_hug_W62.label)
    axes[ai,aj].scatter(qtz_hug_L3.rho,qtz_hug_L3.P,marker='+',color='blue',label=qtz_hug_L3.label)
    axes[ai,aj].plot(qtz_hug_KD.rho,qtz_hug_KD.P,color='green',label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.rho,fs_hug_RTK19.P,'-',color='red',label=fs_hug_RTK19.label)

# silicadb
if 1:
    ind = np.where(qtz_hug_KD.P>50)[0]
    axes[ai,aj].plot(qtz_hug_KD.rho[ind],qtz_hug_KD.P[ind],'g-',linewidth=3,label='KD90 liq fit')
    axes[ai,aj].errorbar(shockdat_list[L83qtzhug['id']].rho,shockdat_list[L83qtzhug['id']].P,fmt='go',mfc='white',\
                 yerr=shockdat_list[L83qtzhug['id']].P_std,label=L83qtzhug['lab'])
    axes[ai,aj].errorbar(shockdat_list[W62qtz['id']].rho,shockdat_list[W62qtz['id']].P,fmt='gx',label=W62qtz['lab'])
    #plt.errorbar(shockdat_list[A23qtz['id']].rho,shockdat_list[A23qtz['id']].P,fmt='gD',mfc='white',label=A23qtz['lab'],\
    #            xerr=shockdat_list[A23qtz['id']].rho_std,yerr=shockdat_list[A23qtz['id']].P_std)
    axes[ai,aj].errorbar(shockdat_list[T20qtz['id']].rho,shockdat_list[T20qtz['id']].P,fmt='<',mfc='white',color='grey',\
                 xerr = shockdat_list[T20qtz['id']].rho_std, yerr=shockdat_list[T20qtz['id']].P_std, label=T20qtz['lab'])

    axes[ai,aj].errorbar(shockdat_list[K13qtz['id']].rho,shockdat_list[K13qtz['id']].P,fmt='gd',markersize=3,\
             xerr=shockdat_list[K13qtz['id']].rho_std,yerr=shockdat_list[K13qtz['id']].P_std,label=K13qtz['lab'])
    ind = np.where(fs_hug_RTK19.P > 74.)[0]
    axes[ai,aj].plot(fs_hug_RTK19.rho[ind],fs_hug_RTK19.P[ind],'b--',label='R19 FS',linewidth=1)
    axes[ai,aj].errorbar(shockdat_list[L03qtz['id']].rho,shockdat_list[L03qtz['id']].P,fmt='r>',label=L03qtz['lab'])

axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.rho,NewEOS.hugo.P,'-.',color='orange',label='ANEOS FS Hugoniot')

axes[ai,aj].set_ylim(0,200)
axes[ai,aj].set_xlim(4.,5.5)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Pressure - Temperature Data vs. ANEOS
ai=0
aj=1

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.T,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.P,NewEOS.hugo.T,'-.',color='orange',label='ANEOS FS Hugoniot')

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Pl,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
    axes[ai,aj].plot(NewEOS.mc.Ps,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

if 0:
    axes[ai,aj].plot(qtz_hug_KD.P,qtz_hug_KD.T/1.e3,color='green',label='Qtz Brygoo et al. 2015')
    axes[ai,aj].plot(fs_hug_Hicks06.P,fs_hug_Hicks06.T/1.e3,'-',color='red',label=fs_hug_Hicks06.label)
    axes[ai,aj].scatter(qtz_hug_LA80.P,qtz_hug_LA80.T/1.E3,marker='+',color='green',label='Qtz Lyzenga & Ahrens 1980')
    axes[ai,aj].plot(stv_mc_MM15.Pl,stv_mc_MM15.T/1.E3,color='cyan',label=stv_mc_MM15.label)

# silicadb
if 1:
    axes[ai,aj].errorbar(shockdat_list[B88qtz['id']].P,shockdat_list[B88qtz['id']].T,fmt='go',mfc='white',yerr=shockdat_list[B88qtz['id']].T_std,label=B88qtz['lab'])
    axes[ai,aj].errorbar(shockdat_list[L04qtz['id']].P,shockdat_list[L04qtz['id']].T,fmt='gv',mfc='white',xerr=shockdat_list[L04qtz['id']].P_std,yerr=shockdat_list[L04qtz['id']].T_std,label=L04qtz['lab'])
    axes[ai,aj].errorbar(shockdat_list[S81qtz['id']].P,shockdat_list[S81qtz['id']].T,fmt='gx',label=S81qtz['lab'])
    # Hicks decaying shock
    axes[ai,aj].plot(qtz_hug_Hicks06.P,qtz_hug_Hicks06.T,'g+:',label='H06 qtz')
    ind = np.where((shockdat_list[B88qtz['id']].P >75) & (shockdat_list[B88qtz['id']].P <110))
    qtztpfit = np.polyfit(shockdat_list[B88qtz['id']].P[ind],shockdat_list[B88qtz['id']].T[ind],2)
    # Hicks decaying shock
    #plt.plot(fs_hug_Hicks06.P,fs_hug_Hicks06.T,'b+:',label='H06 FS')
    axes[ai,aj].errorbar(fs_hug_Hicks06.P,fs_hug_Hicks06.T,fmt='b+:',\
                 yerr=fs_hug_Hicks06.T_std,label='H06 FS')
    print(shockdat_list[B88qtz['id']].P[ind],shockdat_list[B88qtz['id']].T[ind])
    print('quartz Boslough P-T fit 75-110 GPa: ',qtztpfit)
    #plt.plot(shockdat_list[B88qtz['id']].P[ind],qtztpfit[1]+qtztpfit[0]*shockdat_list[B88qtz['id']].P[ind],'g:',label='Likely d-NiAs')
    qtzptmp = np.arange(75,110,1)
    qtzttmp = qtztpfit[2]+qtztpfit[1]*qtzptmp+qtztpfit[0]*qtzptmp*qtzptmp
    axes[ai,aj].plot(qtzptmp,qtzttmp/1.E3,'g--',label='CaCl2 region P-T fit')        
    axes[ai,aj].errorbar(shockdat_list[A23qtz['id']].P,shockdat_list[A23qtz['id']].T,fmt='gD',mfc='white',xerr=shockdat_list[A23qtz['id']].P_std,yerr=shockdat_list[A23qtz['id']].T_std,label=A23qtz['lab'])
    axes[ai,aj].errorbar(shockdat_list[B88fs['id']].P,shockdat_list[B88fs['id']].T,fmt='bo',mfc='white',yerr=shockdat_list[B88fs['id']].T_std,label=B88fs['lab'])

    axes[ai,aj].errorbar(shockdat_list[S80fs['id']].P,shockdat_list[S80fs['id']].T,fmt='bv',xerr=shockdat_list[S80fs['id']].P_std,yerr=shockdat_list[S80fs['id']].T_std,label=S80fs['lab'])
    axes[ai,aj].errorbar(shockdat_list[S81fs['id']].P,shockdat_list[S81fs['id']].T,fmt='bx',label=S81fs['lab'])
    axes[ai,aj].errorbar(shockdat_list[L03qtz['id']].P,shockdat_list[L03qtz['id']].T,fmt='r>',label=L03qtz['lab'])

    
axes[ai,aj].set_xlim(0,200)
axes[ai,aj].set_ylim(0,10000)
axes[ai,aj].set_ylabel('Temperature (1000 K)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Particle Velocity vs. Shock Velocity; Data vs. ANEOS
ai=1
aj=0

axes[ai,aj].plot(NewEOS.hug.up,NewEOS.hug.us,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.up,NewEOS.hugo.us,'-.',color='orange',label='ANEOS FS Hugoniot')

if 0:
    axes[ai,aj].scatter(qtz_hug_W62.up,qtz_hug_W62.us,marker='+',color='purple',label=qtz_hug_W62.label)
    axes[ai,aj].scatter(qtz_hug_L3.up,qtz_hug_L3.us,marker='+',color='blue',label=qtz_hug_L3.label)
    axes[ai,aj].plot(qtz_hug_KD.up,qtz_hug_KD.us,color='green',label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.up,fs_hug_RTK19.us,'-',color='red',label=fs_hug_RTK19.label)

# silicadb
if 1:
    mqtzup = deepcopy(np.asarray(qtzdf['up'][indqtz]))
    mqtzus = deepcopy(np.asarray(qtzdf['us'][indqtz]))
    #inda = np.where((mqtzus>8) & (mqtzus<9.23))[0]
    inda = np.where((mqtzus>6))[0]
    axes[ai,aj].errorbar(mqtzup[inda],mqtzus[inda],fmt='g+',label='M80 qtz included')

    axes[ai,aj].errorbar(shockdat_list[K13qtz['id']].up,shockdat_list[K13qtz['id']].us,fmt='gd',markersize=6,\
                 xerr=shockdat_list[K13qtz['id']].up_std,yerr=shockdat_list[K13qtz['id']].us_std,label=K13qtz['lab'])
    axes[ai,aj].errorbar(shockdat_list[W62fs['id']].up,shockdat_list[W62fs['id']].us,fmt='bx',label=W62fs['lab'])

    axes[ai,aj].errorbar(shockdat_list[L83qtzhug['id']].up,shockdat_list[L83qtzhug['id']].us,fmt='go',mfc='white',\
                 yerr=shockdat_list[L83qtzhug['id']].us_std,label=L83qtzhug['lab'])
    axes[ai,aj].errorbar(shockdat_list[W62qtz['id']].up,shockdat_list[W62qtz['id']].us,fmt='gx',label=W62qtz['lab'])
    #plt.errorbar(shockdat_list[L03qtz['id']].rho,shockdat_list[L03qtz['id']].P,fmt='r>',label=L03qtz['lab'])
    #plt.errorbar(shockdat_list[T20qtz['id']].up,shockdat_list[T20qtz['id']].us,fmt='g<',mfc='white',\
    #             xerr = shockdat_list[T20qtz['id']].up_std, yerr=shockdat_list[T20qtz['id']].us_std, label=T20qtz['lab'])
    axes[ai,aj].errorbar(shockdat_list[T20qtz['id']].up,shockdat_list[T20qtz['id']].us,fmt='<',mfc='white',color='grey',\
                  label=T20qtz['lab'])
    #print(shockdat_list[T20qtz['id']].up,shockdat_list[T20qtz['id']].us)
    #res = np.polyfit(shockdat_list[T20qtz['id']].up,shockdat_list[T20qtz['id']].us,1)
    #print(res)
    axes[ai,aj].errorbar(shockdat_list[L04qtz['id']].up,shockdat_list[L04qtz['id']].us,fmt='gv',mfc='white',\
             xerr=shockdat_list[L04qtz['id']].up_std,yerr=shockdat_list[L04qtz['id']].us_std,label=L04qtz['lab'])
    axes[ai,aj].errorbar(fs_hug_Hicks06.up,fs_hug_Hicks06.us,fmt='b+:',\
                label='H06 FS')
    axes[ai,aj].errorbar(shockdat_list[B19fs['id']].up,shockdat_list[B19fs['id']].us,fmt='b>',mfc='white',label=B19fs['lab'],\
                 xerr=shockdat_list[B19fs['id']].up_std,yerr=shockdat_list[B19fs['id']].us_std)
    axes[ai,aj].errorbar(shockdat_list[L03qtz['id']].up,shockdat_list[L03qtz['id']].us,fmt='r>',label=L03qtz['lab'])

axes[ai,aj].set_xlim(0,6)
axes[ai,aj].set_ylim(0,12)
axes[ai,aj].set_xlabel('Particle Velocity (km/s)',size=sfont)
axes[ai,aj].set_ylabel('Shock Velocity (km/s)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Sp. Entropy vs Pressure; Calculated from experiments vs. ANEOS
ai=1
aj=1

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.S*1.e3,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.P,NewEOS.hugo.S*1.e3,'-.',color='orange',label='ANEOS FS Hugoniot')

if 0:
    # Kraus entropy fit only good to 800 GPa
    ind = np.where(qtz_hug_KD.P <= 800.)
    axes[ai,aj].plot(qtz_hug_KD.P[ind],qtz_hug_KD.S[ind],color='green',label='Qtz Kraus et al. 2012')
    Kraus_serr = (127.3+0.328*qtz_hug_KD.P[ind]-0.000193*np.power(qtz_hug_KD.P[ind],2))/1.e3 # kJ/K/kg
    axes[ai,aj].fill_between(qtz_hug_KD.P[ind],(qtz_hug_KD.S[ind]-Kraus_serr),(qtz_hug_KD.S[ind]+Kraus_serr),color=stylezfit[1],alpha=0.2)

# silicadb
if 1:
    entropyshift = 220.
    penter = 111.
    kqtzparr = np.arange(100,500) # GPa
    kqtzsarr = 1793+22.4*kqtzparr-0.0426*np.power(kqtzparr,2)+4.319e-05*np.power(kqtzparr,3)-1.742e-8*np.power(kqtzparr,4)
    kqtzserrarr = 127.3+0.328*kqtzparr-0.000193*np.power(kqtzparr,2)
    ind = np.where(kqtzparr >= penter)[0]
    axes[ai,aj].plot(kqtzparr[ind],kqtzsarr[ind]/1.e3,'g-',linewidth=1,label='K12 qtz Hugoniot-liq region')
    axes[ai,aj].plot(kqtzparr[ind],(kqtzsarr[ind]-entropyshift)/1.e3,'g-',linewidth=3,label='K12 Qtz Hugoniot shift')
    axes[ai,aj].fill_between(kqtzparr[ind],(kqtzsarr[ind]-entropyshift-kqtzserrarr[ind])/1.e3,(kqtzsarr[ind]-entropyshift+kqtzserrarr[ind])/1.e3,color=stylezfit[1],alpha=0.2)
    axes[ai,aj].errorbar(shockdat_list[L03qtz['id']].P,shockdat_list[L03qtz['id']].S,fmt='r>',label=L03qtz['lab'])
   
    

axes[ai,aj].set_ylim(0,8.)
axes[ai,aj].set_xlim(0,250)
axes[ai,aj].set_ylabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Pressure-Density Data vs. ANEOS   LOW PRESSURE + LIQUID HUGONIOT
ai=2
aj=0
# ANEOS REFERENCE POINT (STP) Hugoniots
axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.rho,NewEOS.hugo.P,'-.',color='orange',label='ANEOS FS Hugoniot')

if 0:
    axes[ai,aj].scatter(qtz_hug_W62.rho,qtz_hug_W62.P,marker='+',color='purple',label=qtz_hug_W62.label)
    axes[ai,aj].scatter(qtz_hug_L3.rho,qtz_hug_L3.P,marker='+',color='blue',label=qtz_hug_L3.label)
    axes[ai,aj].plot(qtz_hug_KD.rho,qtz_hug_KD.P,color='green',label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.rho,fs_hug_RTK19.P,'-',color='red',label=fs_hug_RTK19.label)

# silicadb
if 1:
    axes[ai,aj].errorbar(shockdat_list[L83qtzhug['id']].rho,shockdat_list[L83qtzhug['id']].P,fmt='go',mfc='white',\
                 yerr=shockdat_list[L83qtzhug['id']].P_std,label=L83qtzhug['lab'])
    axes[ai,aj].errorbar(shockdat_list[W62qtz['id']].rho,shockdat_list[W62qtz['id']].P,fmt='gx',label=W62qtz['lab'])
    #plt.errorbar(shockdat_list[A23qtz['id']].rho,shockdat_list[A23qtz['id']].P,fmt='gD',mfc='white',label=A23qtz['lab'],\
    #            xerr=shockdat_list[A23qtz['id']].rho_std,yerr=shockdat_list[A23qtz['id']].P_std)
    axes[ai,aj].errorbar(shockdat_list[T20qtz['id']].rho,shockdat_list[T20qtz['id']].P,fmt='<',mfc='white',color='grey',\
                 xerr = shockdat_list[T20qtz['id']].rho_std, yerr=shockdat_list[T20qtz['id']].P_std, label=T20qtz['lab'])
    #print(shockdat_list[T20qtz['id']].up,shockdat_list[T20qtz['id']].us)
    #res = np.polyfit(shockdat_list[T20qtz['id']].up,shockdat_list[T20qtz['id']].us,1)
    #print(res)
    axes[ai,aj].errorbar(shockdat_list[L04qtz['id']].rho,shockdat_list[L04qtz['id']].P,fmt='gv',mfc='white',\
             xerr=shockdat_list[L04qtz['id']].rho_std,yerr=shockdat_list[L04qtz['id']].P_std,label=L04qtz['lab'])

    ind = np.where(qtz_hug_KD.P>112)[0]
    axes[ai,aj].plot(qtz_hug_KD.rho[ind],qtz_hug_KD.P[ind],'g-',linewidth=3,label='KD90 liq fit')
    axes[ai,aj].errorbar(qtzdf['rho'][indqtz],qtzdf['P'][indqtz],fmt='+',color='grey',label='M80 qtz')

    
axes[ai,aj].set_ylim(-5,200)
axes[ai,aj].set_xlim(2.5,6.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Pressure-Density Data vs. ANEOS   EXTRAPOLATED HIGH PRESSURE
ai=2
aj=1
axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.rho,NewEOS.hugo.P,'-.',color='orange',label='ANEOS FS Hugoniot')

if 0:
    axes[ai,aj].plot(qtz_hug_KD.rho,qtz_hug_KD.P,color='green',label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.rho,fs_hug_RTK19.P,'-',color='red',label=fs_hug_RTK19.label)

# silicadb
if 1:
    axes[ai,aj].errorbar(shockdat_list[W62fs['id']].rho,shockdat_list[W62fs['id']].P,fmt='bx',label=W62fs['lab'])
    axes[ai,aj].errorbar(shockdat_list[L83fshug['id']].rho,shockdat_list[L83fshug['id']].P,fmt='bo',label=L83fshug['lab'])
    axes[ai,aj].errorbar(shockdat_list[B19fs['id']].rho,shockdat_list[B19fs['id']].P,fmt='b>',xerr=shockdat_list[B19fs['id']].rho_std,\
                 yerr=shockdat_list[B19fs['id']].P_std,mfc='white',label=B19fs['lab'])
    axes[ai,aj].errorbar(shockdat_list[T18fs['id']].rho,shockdat_list[T18fs['id']].P,fmt='b<',mfc='white',\
                 xerr=shockdat_list[T18fs['id']].rho_std,yerr=shockdat_list[T18fs['id']].P_std,label=T18fs['lab'])

    #plt.errorbar(np.asarray(shockdat_list[A23fs['id']].rho)/1000.,shockdat_list[A23fs['id']].P,fmt='bD',mfc='white',label=A23fs['lab'])
    #plt.errorbar(np.asarray(shockdat_list[A23fs['id']].rho),shockdat_list[A23fs['id']].P,fmt='bD',mfc='white',label=A23fs['lab'],\
    #            xerr=shockdat_list[A23fs['id']].rho_std,yerr=shockdat_list[A23fs['id']].P_std)
    axes[ai,aj].errorbar(fsdf['rho'][indfs],fsdf['P'][indfs],fmt='+',color='grey',label='M80 FS excluded')
    indhigh = np.where(fsdf['P'][indfs] > 74)[0]
    axes[ai,aj].errorbar(fsdf['rho'][indfs][indhigh],fsdf['P'][indfs][indhigh],fmt='b+',label='M80 FS included')

    #axes[ai,aj].plot(fitfslrarr,fitfslparr,'b-',linewidth=2,label='cubic fit liquid')
    #axes[ai,aj].plot(fitfssrarr,fitfssparr,'b-',label='linear fit solid')
    ind = np.where(fs_hug_RTK19.P > 74.)[0]
    axes[ai,aj].plot(fs_hug_RTK19.rho[ind],fs_hug_RTK19.P[ind],'b--',label='R19 FS',linewidth=1)


axes[ai,aj].set_ylim(0,300)
axes[ai,aj].set_xlim(3.,6.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

plt.savefig(MODELNAME+'-Hugoniots.pdf')

Requested ANEOS Hugoniot initial state from (THUG,RHUG) (red curves in plots below):
Initial Pressure (GPa): 0.0001
Initial Temperature (K): 298.0
Initial Sp. Entropy (MJ/K/kg): 0.00017886
Initial Sp. Internal Energy (MJ/kg): 0.043369
Initial Density (g/cm3): 2.2
[ 75.9  85.9  92.5  99.3 107.8 109.7] [4440. 4850. 5430. 5630. 5900. 6100.]
quartz Boslough P-T fit 75-110 GPa:  [-3.39692149e-01  1.11592149e+02 -2.10423906e+03]

print(shockdat_list[L03qtz['id']].S)

[3.365 3.527 3.621 3.712 3.817 3.839]

# Display the initial conditions for the Hugoniot in ANEOS.OUTPUT    
print('\nRequested ANEOS Hugoniot initial state from (THUG,RHUG) (red curves in plots below):')
print('Initial Pressure (GPa):',NewEOS.anhug.P[0])
print('Initial Temperature (K):',NewEOS.anhug.T[0])
print('Initial Sp. Entropy (MJ/K/kg):',NewEOS.anhug.S[0])
print('Initial Sp. Internal Energy (MJ/kg):',NewEOS.anhug.U[0])
print('Initial Density (g/cm3):',NewEOS.anhug.rho[0])
#

#========== PLOTTING STYLES FOR LOTS OF PLOTS =================
styleqmd = ['^','red',5,'Fo QMD calc (Root et al. 2018)']
stylezdat = ['^','black',5,'Fo Z data (Root et al. 2018)']
stylegg = ['+','green',5,'Fo Gas gun data\n(Mosenfelder et al. 2007)']
styleliq = ['^','orange',7,'Liquid Fo data corrected\n(Thomas & Asimow 2013)']
styleliqhug = ['-','orange',1,'ANEOS Liquid 2273K Hugoniot']
stylezfit = ['-','green',1,'Fo Z fit (Davies et al. 2019)']
styleanhug = ['-','red',1,'ANEOS Glass-tension Hugoniot']
stylestphug = ['-','blue',1,'Tabulated STP Hugoniot']
styleanphase = ['-','black',1,'ANEOS Phase Boundaries']
stylephasept = ['o','black',6]



#======================LOTS OF PLOTS============================
# plot the Hugoniots for comparison
sfont = 15
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(16,12))
plt.subplots_adjust(wspace=0.4)

#-------------
# Pressure-Density Data vs. ANEOS
ai=0
aj=0
axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.rho,NewEOS.hugo.P,'-.',color='orange',label='ANEOS FS Hugoniot')

if 0:
    axes[ai,aj].plot(qtz_hug_KD.rho,qtz_hug_KD.P,color='green',label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.rho,fs_hug_RTK19.P,'-',color='red',label=fs_hug_RTK19.label)


axes[ai,aj].set_ylim(0,1400)
axes[ai,aj].set_xlim(3.,9.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Pressure - Temperature Data vs. ANEOS
ai=0
aj=1

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.T/1.E3,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label='ANEOS STP Hugoniot')
axes[ai,aj].plot(NewEOS.hugo.P,NewEOS.hugo.T/1.E3,'-.',color='lime',label='ANEOS Glass Hugoniot')

if 0:
    axes[ai,aj].plot(qtz_hug_KD.P,qtz_hug_KD.T/1.e3,color='green',label='Qtz Brygoo et al. 2015')
    axes[ai,aj].plot(fs_hug_Hicks06.P,fs_hug_Hicks06.T/1.e3,'-',color='red',label=fs_hug_Hicks06.label)

    axes[ai,aj].scatter(qtz_hug_LA80.P,qtz_hug_LA80.T/1.E3,marker='+',color='green',label='Qtz Lyzenga & Ahrens 1980')
    axes[ai,aj].plot(stv_mc_MM15.Pl,stv_mc_MM15.T/1.E3,color='cyan',label=stv_mc_MM15.label)


axes[ai,aj].set_xlim(0,1200)
axes[ai,aj].set_ylim(0,70)
axes[ai,aj].set_ylabel('Temperature (1000 K)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)


#-------
# pressure - temperature melt curve
ai=1
aj=0
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Pl,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
    axes[ai,aj].plot(NewEOS.mc.Ps,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.T,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.P[0],NewEOS.hug.T[0],'o',color=stylestphug[1])

axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(0,15000)
axes[ai,aj].set_xlim(-5,600)

#axes[ai,aj].plot(peridotiteliquidusdatarr.iloc[:,0].values,peridotiteliquidusdatarr.iloc[:,1].values,linestyle='-',color='orange',label='Peridotite Liquidus')
#axes[ai,aj].plot(peridotitesolidusdatarr.iloc[:,0].values,peridotitesolidusdatarr.iloc[:,1].values,linestyle=':',color='orange',label='Peridotite Solidus\n(Fiquet et al. 2010)')

if 1:
    #axes[ai,aj].plot(stv_mc_MM15.Pl,stv_mc_MM15.T,color='cyan',label=stv_mc_MM15.label)
    axes[ai,aj].plot(qtz_mc_BR98.Pl,qtz_mc_BR98.T,color='pink',label=qtz_mc_BR98.label)

if 1:
    wcol = "#653700"
    wcol2 = "#FBDD7E"
    wcol3 = "#F97306"
    wop = 0.5
    wop2=0.75
    wop3=0.25
    ind = np.arange(0,11)
    print(stsmc.P[ind])
    axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],label='melt curve exp',alpha=wop,color=wcol)
    ind = np.arange(12,16)
    print(stsmc.P[ind])
    axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],alpha=wop,color=wcol)
    ind = np.arange(17,21)
    print(stsmc.P[ind])
    axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],alpha=wop,color=wcol)
    ind = np.arange(22,24)
    print(stsmc.P[ind])
    axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],alpha=wop2,color=wcol2,label='melt curve calc')
    ind = np.arange(25,28)
    print(stsmc.P[ind])
    axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],alpha=wop2,color=wcol2)
    # Millot stv melt curve crossing
    axes[ai,aj].fill_between([470.,530.],[8000.,8000.],[8600.,8600.],alpha=0.5,color='blue',label='MM15')
    

if 0:
    ind = np.where(stsmcsmooth.P <=Pt3)[0]
    axes[ai,aj].plot(stsmcsmooth.P[ind],stsmcsmooth.T[ind],'m-',linewidth=3,label='Proposed melt curve')
    ind = np.where(stsmcsmooth.P >=Pt3)[0]
    axes[ai,aj].plot(stsmcsmooth.P[ind],stsmcsmooth.T[ind],'m:',linewidth=3)#,label='Proposed melt curve')



axes[ai,aj].legend()

#-------
ai=1
aj=1
#axes[ai,aj].axis('off')

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.S*1.e3,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.P,NewEOS.hugo.S*1.e3,'-.',color='orange',label='ANEOS FS Hugoniot')

if 0:
    # Kraus entropy fit only good to 800 GPa
    ind = np.where(qtz_hug_KD.P <= 800.)
    axes[ai,aj].plot(qtz_hug_KD.P[ind],qtz_hug_KD.S[ind],color='green',label='Qtz Kraus et al. 2012')
    Kraus_serr = (127.3+0.328*qtz_hug_KD.P[ind]-0.000193*np.power(qtz_hug_KD.P[ind],2))/1.e3 # kJ/K/kg
    axes[ai,aj].fill_between(qtz_hug_KD.P[ind],(qtz_hug_KD.S[ind]-Kraus_serr),(qtz_hug_KD.S[ind]+Kraus_serr),color=stylezfit[1],alpha=0.2)

# silicadb
if 1:
    entropyshift = 220.
    penter = 111.
    kqtzparr = np.arange(100,500) # GPa
    kqtzsarr = 1793+22.4*kqtzparr-0.0426*np.power(kqtzparr,2)+4.319e-05*np.power(kqtzparr,3)-1.742e-8*np.power(kqtzparr,4)
    kqtzserrarr = 127.3+0.328*kqtzparr-0.000193*np.power(kqtzparr,2)
    ind = np.where(kqtzparr >= penter)[0]
    axes[ai,aj].plot(kqtzparr[ind],kqtzsarr[ind]/1.e3,'g-',linewidth=1,label='K12 qtz Hugoniot-liq region')
    axes[ai,aj].plot(kqtzparr[ind],(kqtzsarr[ind]-entropyshift)/1.e3,'g-',linewidth=3,label='K12 Qtz Hugoniot shift')
    axes[ai,aj].fill_between(kqtzparr[ind],(kqtzsarr[ind]-entropyshift-kqtzserrarr[ind])/1.e3,(kqtzsarr[ind]-entropyshift+kqtzserrarr[ind])/1.e3,color=stylezfit[1],alpha=0.2)
   
    

axes[ai,aj].set_ylim(0,5.)
axes[ai,aj].set_xlim(0,200)
axes[ai,aj].set_ylabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

plt.savefig(MODELNAME+'-3panel.pdf')

Requested ANEOS Hugoniot initial state from (THUG,RHUG) (red curves in plots below):
Initial Pressure (GPa): 0.0001
Initial Temperature (K): 298.0
Initial Sp. Entropy (MJ/K/kg): 0.00017886
Initial Sp. Internal Energy (MJ/kg): 0.043369
Initial Density (g/cm3): 2.2
[13.7  14.5  20.   23.75 28.   34.   36.   37.5  40.   42.   46.  ]
[55. 60. 66. 74.]
[100. 102. 114. 117.]
[135. 155.]
[217.5 218.5   nan]

print('ANEOS Calculated Critical point:')
print('   Temperature (K):       ',NewEOS.cp.T)
print('   Pressure (GPa):        ', NewEOS.cp.P)
print('   Density (g/cm^3):      ',NewEOS.cp.rho)
print('   Sp. Energy (MJ/kg):    ', NewEOS.cp.U)
print('   Sp. Entropy (kJ/K/kg): ', NewEOS.cp.S*1.e3)
print('\n')
print('ANEOS Calculated Values at 1 bar:')
print('   Temperature at melting (K):                     ', NewEOS.onebar.Tmelt)
print('   Liquid density at melting (g/cm^3):             ', NewEOS.onebar.rhocm)
print('   Sp. Entropy at inc. melting Sim (kJ/K/kg):      ', NewEOS.onebar.Sim*1.E3)
print('   Sp. Entropy at com. melting Scm (kJ/K/kg):      ', NewEOS.onebar.Scm*1.E3)
print('   Temperature at boiling (K):                     ', NewEOS.onebar.Tvap)
print('   Liquid density at boiling (g/cm^3):             ', NewEOS.onebar.rhoiv)
print('   Sp. Entropy at inc. vaporization Siv (kJ/K/kg): ', NewEOS.onebar.Siv*1.E3)
print('   Sp. Entropy at com. vaporization Scv (kJ/K/kg): ', NewEOS.onebar.Scv*1.E3)
print('\n')
if NewEOS.tp.T > 0:
    print('ANEOS Calculated Triple Point:')
    print('   Temperature (K):                        ',NewEOS.tp.T)
    print('   Pressure (GPa):                         ', NewEOS.tp.P)
    print('   Solid density (g/cm^3):                 ',NewEOS.tp.rhos)
    print('   Liquid density (g/cm^3):                ',NewEOS.tp.rhol)
    print('   Vapor density (g/cm^3):                 ',NewEOS.tp.rhov)
    print('   Sp. Entropy at inc. melt Sim (kJ/K/kg): ', NewEOS.tp.Sim*1.E3)
    print('   Sp. Entropy at com. melt Scm (kJ/K/kg): ', NewEOS.tp.Scm*1.E3)
    print('   Sp. Entropy at inc. vap. Siv (kJ/K/kg): ', NewEOS.tp.Siv*1.E3)
    print('   Sp. Entropy at com. vap. Scv (kJ/K/kg): ', NewEOS.tp.Scv*1.E3)
else:
    print('ANEOS: No melt curve or triple point in this calculation.')
print('\n')
print('Starting at reference conditions (rho,T): ',NewEOS.R0REF,NewEOS.T0REF, ' (g/cm^3, K)')
print('ANEOS model conditions for shock-induced vaporization on release to 1 bar:')
print('                     :         IM        CM        IV      50%V        CV')
print('Shock Pressure (GPa) : {:10.2f}{:10.2f}{:10.2f}{:10.2f}{:10.2f}'.format(
    np.interp(NewEOS.onebar.Sim,NewEOS.hug.S,NewEOS.hug.P),
    np.interp(NewEOS.onebar.Scm,NewEOS.hug.S,NewEOS.hug.P),
    np.interp(NewEOS.onebar.Siv,NewEOS.hug.S,NewEOS.hug.P),
    np.interp(NewEOS.onebar.Siv+0.5*(NewEOS.onebar.Scv-NewEOS.onebar.Siv),NewEOS.hug.S,NewEOS.hug.P),
    np.interp(NewEOS.onebar.Scv,NewEOS.hug.S,NewEOS.hug.P)))
print('Imp. Match Vel (km/s): {:10.2f}{:10.2f}{:10.2f}{:10.2f}{:10.2f}'.format(
    np.interp(NewEOS.onebar.Sim,NewEOS.hug.S,NewEOS.hug.up)*2,
    np.interp(NewEOS.onebar.Scm,NewEOS.hug.S,NewEOS.hug.up)*2,
    np.interp(NewEOS.onebar.Siv,NewEOS.hug.S,NewEOS.hug.up)*2,
    np.interp(NewEOS.onebar.Siv+0.5*(NewEOS.onebar.Scv-NewEOS.onebar.Siv),NewEOS.hug.S,NewEOS.hug.up)*2,
    np.interp(NewEOS.onebar.Scv,NewEOS.hug.S,NewEOS.hug.up)*2))
print('ANEOS model conditions for shock-induced vaporization on release to the triple point pressure:')
print('                     :         IM        CM        IV      50%V        CV')
print('Shock Pressure (GPa) : {:10.2f}{:10.2f}{:10.2f}{:10.2f}{:10.2f}'.format(
    np.interp(NewEOS.tp.Sim,NewEOS.hug.S,NewEOS.hug.P),
    np.interp(NewEOS.tp.Scm,NewEOS.hug.S,NewEOS.hug.P),
    np.interp(NewEOS.tp.Siv,NewEOS.hug.S,NewEOS.hug.P),
    np.interp(NewEOS.tp.Siv+0.5*(NewEOS.tp.Scv-NewEOS.tp.Siv),NewEOS.hug.S,NewEOS.hug.P),
    np.interp(NewEOS.tp.Scv,NewEOS.hug.S,NewEOS.hug.P)))
print('Imp. Match Vel (km/s): {:10.2f}{:10.2f}{:10.2f}{:10.2f}{:10.2f}'.format(
    np.interp(NewEOS.tp.Sim,NewEOS.hug.S,NewEOS.hug.up)*2,
    np.interp(NewEOS.tp.Scm,NewEOS.hug.S,NewEOS.hug.up)*2,
    np.interp(NewEOS.tp.Siv,NewEOS.hug.S,NewEOS.hug.up)*2,
    np.interp(NewEOS.tp.Siv+0.5*(NewEOS.tp.Scv-NewEOS.tp.Siv),NewEOS.hug.S,NewEOS.hug.up)*2,
    np.interp(NewEOS.tp.Scv,NewEOS.hug.S,NewEOS.hug.up)*2))

ANEOS Calculated Critical point:
   Temperature (K):        7512.49
   Pressure (GPa):         0.390284
   Density (g/cm^3):       0.507337
   Sp. Energy (MJ/kg):     16.8021
   Sp. Entropy (kJ/K/kg):  5.76996


ANEOS Calculated Values at 1 bar:
   Temperature at melting (K):                      1996.010972480397
   Liquid density at melting (g/cm^3):              1.8810072665433093
   Sp. Entropy at inc. melting Sim (kJ/K/kg):       2.4828055952383483
   Sp. Entropy at com. melting Scm (kJ/K/kg):       2.8135042388169302
   Temperature at boiling (K):                      3212.86
   Liquid density at boiling (g/cm^3):              1.6832
   Sp. Entropy at inc. vaporization Siv (kJ/K/kg):  3.58174
   Sp. Entropy at com. vaporization Scv (kJ/K/kg):  6.80183


ANEOS Calculated Triple Point:
   Temperature (K):                         1996.0
   Pressure (GPa):                          2.90104e-09
   Solid density (g/cm^3):                  2.0225999999999997
   Liquid density (g/cm^3):                 1.881
   Vapor density (g/cm^3):                  7.5663e-09
   Sp. Entropy at inc. melt Sim (kJ/K/kg):  2.4828
   Sp. Entropy at com. melt Scm (kJ/K/kg):  2.8135
   Sp. Entropy at inc. vap. Siv (kJ/K/kg):  2.81349
   Sp. Entropy at com. vap. Scv (kJ/K/kg):  8.45833


Starting at reference conditions (rho,T):  2.2 298.0  (g/cm^3, K)
ANEOS model conditions for shock-induced vaporization on release to 1 bar:
                     :         IM        CM        IV      50%V        CV
Shock Pressure (GPa) :      39.70     50.46     84.62    243.49    725.55
Imp. Match Vel (km/s):       6.05      6.86      9.01     15.92     28.78
ANEOS model conditions for shock-induced vaporization on release to the triple point pressure:
                     :         IM        CM        IV      50%V        CV
Shock Pressure (GPa) :      39.70     50.46     50.46    330.10   1872.74
Imp. Match Vel (km/s):       6.05      6.86      6.86     18.79     47.50

# PHASE DIAGRAM PLOTS
#======================LOTS OF PLOTS============================
# Plot aneos phase boundaries with focus on the vapor curve
sfont = 15
fig, axes = plt.subplots(nrows=4, ncols=2, figsize=(14,24))
plt.subplots_adjust(wspace=0.4)

#------------------------------
# entropy-temperature
ai=0
aj=0
axes[ai,aj].plot(NewEOS.vc.Sl*1.E3,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
axes[ai,aj].plot(NewEOS.vc.Sv*1.E3,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
axes[ai,aj].plot(NewEOS.cp.S*1.E3,NewEOS.cp.T,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Sl*1.E3,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
    axes[ai,aj].plot(NewEOS.mc.Ss*1.E3,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.S*1.E3,NewEOS.hug.T,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.S[0]*1.E3,NewEOS.hug.T[0],'o',color=stylestphug[1])

axes[ai,aj].set_ylim(0,1.5*NewEOS.cp.T)
axes[ai,aj].set_xlim(0,14)
axes[ai,aj].set_xlabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_ylabel('Temperature [K]',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-------------------------
# density-pressure
ai=0
aj=1
axes[ai,aj].loglog(NewEOS.vc.rl,NewEOS.vc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.Pv,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.P,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.rho[0],NewEOS.hug.P[0],'o',color=stylestphug[1])



axes[ai,aj].set_ylim(1.e-5,300)
axes[ai,aj].set_xlim(1.e-1,8)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#------
#entropy-pressure
ai=1
aj=0
axes[ai,aj].semilogy(NewEOS.vc.Sl*1.E3,NewEOS.vc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
axes[ai,aj].plot(NewEOS.vc.Sv*1.E3,NewEOS.vc.Pv,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
axes[ai,aj].plot(NewEOS.cp.S*1.E3,NewEOS.cp.P,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Sl*1.E3,NewEOS.mc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
    axes[ai,aj].plot(NewEOS.mc.Ss*1.E3,NewEOS.mc.Ps,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.S*1.E3,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.S[0]*1.E3,NewEOS.hug.P[0],'o',color=stylestphug[1])

axes[ai,aj].set_xlabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(1.e-10,10000)
axes[ai,aj].set_xlim(0,14)

#-----------
# specific energy - pressure
ai=1
aj=1
axes[ai,aj].semilogy(NewEOS.vc.Ul,NewEOS.vc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
axes[ai,aj].plot(NewEOS.vc.Uv,NewEOS.vc.Pv,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
axes[ai,aj].plot(NewEOS.cp.U,NewEOS.cp.P,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Ul,NewEOS.mc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
    axes[ai,aj].plot(NewEOS.mc.Us,NewEOS.mc.Ps,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.U,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.U[0],NewEOS.hug.P[0],'o',color=stylestphug[1])

axes[ai,aj].set_xlabel('Specific Energy (MJ/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(1.e-11,10000)
axes[ai,aj].set_xlim(-3,NewEOS.cp.U*2.)

#-------
# temperature - pressure
ai=2
aj=0
axes[ai,aj].semilogy(NewEOS.vc.T,NewEOS.vc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
axes[ai,aj].plot(NewEOS.vc.T,NewEOS.vc.Pv,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
axes[ai,aj].plot(NewEOS.cp.T,NewEOS.cp.P,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.T,NewEOS.mc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
    axes[ai,aj].plot(NewEOS.mc.T,NewEOS.mc.Ps,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.T,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.T[0],NewEOS.hug.P[0],'o',color=stylestphug[1])

axes[ai,aj].set_xlabel('Temperature (K)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_xlim(0,9000)
axes[ai,aj].set_ylim(1.e-16,100)

#-------
# pressure - temperature melt curve
ai=2
aj=1
#axes[ai,aj].plot(NewEOS.vc.Pl,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
#axes[ai,aj].plot(NewEOS.vc.Pv,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
#axes[ai,aj].plot(NewEOS.cp.P,NewEOS.cp.T,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Pl,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
    axes[ai,aj].plot(NewEOS.mc.Ps,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.T,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.P[0],NewEOS.hug.T[0],'o',color=stylestphug[1])
axes[ai,aj].plot(NewEOS.hugo.P,NewEOS.hugo.T,'-.',color='orange',label='ANEOS FS Hugoniot')

axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(0,9000)
axes[ai,aj].set_xlim(-5,600)


#-------
# density - pressure melt curve
ai=3
aj=0
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.rho[0],NewEOS.hug.P[0],'o',color=stylestphug[1])

axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(0,300)
axes[ai,aj].set_xlim(2.5,6.5)


#-------
# density -- temperature
ai=3
aj=1
axes[ai,aj].semilogx(NewEOS.vc.rl,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.T,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.T,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
#axes[ai,aj].plot(NewEOS.hug.T[0],NewEOS.hug.P[0],'o',color=stylestphug[1])

axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(0,9000)
axes[ai,aj].set_xlim(0.1,10)



    
#===================================================================
# PLOT EXPERIMENTAL CONSTRAINTS ON THE ANEOS PHASE BOUNDARIES
PLOTDATA=1
datms = 8
if PLOTDATA:
    # entropy-temperature [0,0]
    ai=0
    aj=0
    
    # Connolly JGR 2016 summary of critical points
    axes[ai,aj].add_patch(Rectangle((4.8,5071), 0.5, 1000,  alpha=0.7, facecolor='lime', linewidth=0))
    axes[ai,aj].scatter(1900, -1.e-9, alpha=0.7, marker='s',color='lime', label='Connolly2016 CP')
    
    axes[ai,aj].plot(5.149,5130,'^',markersize=datms,color='black',label='Qtz Kraus et al. 2012')
    # Triple point from Nagahara et al. 1994
    #Triple point of forsterite 5.2e-5 bar, 1890 C (2163 K)
    axes[ai,aj].plot(3.470,2163,'X',color='orange',label='Fo TP Nagahara et al. 1994')

    axes[ai,aj].add_patch(Rectangle((1850, 1.e-10), 600, 3.e-7, alpha=0.2, facecolor='brown', linewidth=0))
    axes[ai,aj].plot([3.507,7.267],[3127,3127],'d',color='red',label='JANAF bp')
    axes[ai,aj].plot([3.507],[1996],'s',color='red',label='MK88 tp')

    
    if 1:
        #plot Connolly calcs
        b = np.argsort(C16highs)
        axes[ai,aj].plot(C16highs[b]/1.e3,C16hight[b],'-',color='black',label='C16 HSC-vWRK vapor curve')
        axes[ai,aj].plot(5.382,6074,'o',color='black',markersize=5)
        axes[ai,aj].plot(C16lows/1.e3,C16lowt,'--',color='black',label='C16 HSC-CSRK vapor curve')
        axes[ai,aj].plot(5.322,5548,'o',color='black',markersize=5)
        #plt.plot(C16sils,C16silt,'-',color='orangered',label='C16 SiO(l)-vWRK')
        #plt.plot(5403,5675,'o',color='orangered')
        # plot boiling point at 1 bar
        # note Kraus S_vap too small
        axes[ai,aj].plot([3.554,7.460],[3130,3130],'kD',label='C16 1-bar boiling pt.')
      
    # density-pressure [0,1]
    ai=0
    aj=1
    # Connolly JGR 2016 summary of critical points
    axes[ai,aj].add_patch(Rectangle((0.5,.13), 0.5, .8,  alpha=0.7, facecolor='lime', linewidth=0))
    axes[ai,aj].scatter(1900, -1.e-9, alpha=0.7, marker='s',color='lime', label='Connolly2016 CP')

    axes[ai,aj].plot(0.4822,0.1414,'+',label='En Xiao & Stixrude 2018',markersize=datms,color='brown')
    axes[ai,aj].plot(0.508,0.132,'^',label='Qtz Kraus et al. 2012',markersize=datms,color='black')
    axes[ai,aj].plot(0.51,0.114,'*',label='Fo Townsend et al. 2020',markersize=datms,color='magenta')

    axes[ai,aj].plot(60.083/27.30,1.e-4,'s',color='red',label='Bourova Richet 1998 qtz')
    axes[ai,aj].plot(60.083/23.45,1.e-4,'d',color='orange',label='Bourova Richet 1998 liq')
    axes[ai,aj].plot(60.083/23.51,1.e-4,'d',color='orange',label='')

    # Caracas Pyrolite DFT
    #axes[ai,aj].add_patch(Rectangle((.3, .08), 0.7, 0.05, alpha=0.5, facecolor='cyan', linewidth=0))
    
    
    # entropy -- pressure [1,0]
    ai=1
    aj=0
    # 0.6689 kJ/kg STP entropy from Robie et al. 1982
    axes[ai,aj].plot(0.669,1.e-4,'+',color='orange',label='STP')
    # 1 bar melting point entropy, enthalpy and temperature from Richet et al. 1993
    axes[ai,aj].plot(3.470,1.e-4,'d',color='orange',label='MP 1 bar Richet et al. 1993')
    # Triple point from Nagahara et al. 1994
    #Triple point of forsterite 5.2e-5 bar, 1890 C (2163 K)
    axes[ai,aj].plot(3.470,5.2e-9,'X',color='orange',label='TP Nagahara et al. 1994')

    axes[ai,aj].plot([3.507,7.267],[1.e-4,1.e-4],'d',color='red',label='JANAF bp')

    
    # temperature -- pressure [2,0]
    ai=2
    aj=0
    
    # Petaev BSE triple system region
    #axes[ai,aj].add_patch(Rectangle((1850, 1.e-10), 600, 3.e-7, alpha=0.2, facecolor='brown', linewidth=0))
    #axes[ai,aj].scatter(1900, -1.e-9, alpha=0.2, marker='s',color='brown', label='Petaev BSE TP')
    
    # Triple point from Nagahara et al. 1994
    #Triple point of forsterite 5.2e-5 bar, 1890 C (2163 K)
    axes[ai,aj].plot(2163.,5.2e-9,'X',color='orange',label='Fo TP Nagahara et al. 1994')
    # Melt curve points
    # axes[ai,aj].plot(2174,1.e-4,'d',color='orange',label='Richet et al. 1993')
    # axes[ai,aj].plot(2550,13.,'v',color='orange',label='de Koker et al. 2009')
    # Nagahara et al. 1994 vapor curve of forsterite
    tvapfo = np.arange(46)*10.+1400.+273. # in C to K
    #pvapfo = np.exp(-77955.0/tvapfo+27.7)
    pvapfonagaqhara = np.exp(169.-543.e3/(8.314*(tvapfo))) 
    #print(pvapfonagahara)
    #print(tvapfo)
    #print(pvapfo)
    #axes[3,0].plot(tvapfo+273.,pvapfo*1.e5/1.e9,linestyle='-',color='orange',label='Sublimation curve\nMysen & Kushiro 1988')
    #axes[3,0].plot(tvapfo,pvapfonagahara*1.e5/1.e9,linestyle='-',color='green',label='Sublimation curve\nNagahara et al. 1994')
    #print(pvapfonagahara)
    #tt=np.zeros(7)
    #tt.fill(273.)
    #tt=tt+np.asarray([1400.,1500.,1600.,1700.,1725.,1760.,1800.])
    #axes[ai,aj].plot(tt,np.asarray([8.79E-9,5.98E-8,3.85E-7,2.7E-6,4.21E-6,9.58E-6,1.52E-5])/1.e4,'-x',color='orange',label='Fo Nagahara et al. 1994')

    #Centolanzi and Chapman Tektite vapor curve
    tvaptek = np.arange(1400)+1700. # in K
    pvaptek = np.exp(-57400.0/tvaptek+18.5)*101325. # in atm to Pa
    axes[ai,aj].plot(tvaptek,pvaptek/1.e9,linestyle='-',color='cyan',label='Vapor curve tektites\nCentolanzi & Chapman 1966')

    # Misha Petaev BSE vapor curve 50% condensed
    #axes[ai,aj].scatter([2500., 3000., 3500., 4000., 4500.],[0.5e-6,2.e-5,3.e-4,20e-4,100.e-4],marker='+',label='Petaev BSE 50% condensation')
    #axes[ai,aj].plot(6110,0.114,'*',label='Fo Townsend et al. 2020',markersize=datms,color='magenta')

    # Mysen and Kushiro 1988 SiO2 vapor curve and TPs
    # not sure this is pure SiO2 - may be mixed with H2 - check
    #axes[ai,aj].add_patch(Rectangle((1608+273, 4.e-11), 93, 1.3e-10, alpha=0.7, facecolor='green', linewidth=0))
    #axes[ai,aj].scatter(1900, -1.e-9, alpha=0.7, marker='s',color='green', label='Mysen1988 TP')
    # fitted temps between 2000 and 1880 - this is SiO2+H2
    #MKvap_t = np.arange(1880,2000,20) # K Fig 5
    #MKvap_p = np.exp(-50011/MKvap_t + 14.2)*1.e5/1.e9 # GPa
    #axes[ai,aj].plot(MKvap_t,MKvap_p,color='red',label='Mysen1988 SiO2+H2 VC')
    #print(MKvap_t,MKvap_p)
    # Figure 2 pure SiO2 Hidalgo 1960 equation 1660-1920 K
    MKvap_t = np.arange(1660,1920,20) # K Fig 2
    MKvap_p = np.exp(-57800/MKvap_t + 18.4)*1.e5/1.e9 # GPa
    axes[ai,aj].plot(MKvap_t,MKvap_p,color='red',label='MK88/H60 SiO2 VC')
    
    # Connolly JGR 2016 summary of critical points
    axes[ai,aj].add_patch(Rectangle((5071, .13), 1000, .8, alpha=0.7, facecolor='lime', linewidth=0))
    axes[ai,aj].scatter(1900, -1.e-9, alpha=0.7, marker='s',color='lime', label='Connolly2016 CP')
    
    # pressure- temperature [2,1]
    ai=2
    aj=1
    # Melt curve points
    axes[ai,aj].plot(stv_mc_MM15.Pl,stv_mc_MM15.T,label=stv_mc_MM15.label)
    axes[ai,aj].plot(qtz_mc_BR98.Pl,qtz_mc_BR98.T,label=qtz_mc_BR98.label)

    #axes[ai,aj].scatter(qtz_hug_LA80.P,qtz_hug_LA80.T,label=qtz_hug_LA80.label)
    #axes[ai,aj].plot(qtz_hug_KD.P,qtz_hug_KD.T,color='green',label='Qtz Brygoo et al. 2015')


    #axes[ai,aj].plot(Mosenfelderdatarr.iloc[:,0].values,Mosenfelderdatarr.iloc[:,1].values,linestyle='-',color='brown',linewidth=1,label='Forsterite Melt Curve\n(Mosenfelder et al. 2007)')
    #axes[ai,aj].plot(peridotiteliquidusdatarr.iloc[:,0].values,peridotiteliquidusdatarr.iloc[:,1].values,linestyle='-',color='orange',label='Peridotite Liquidus')
    #axes[ai,aj].plot(peridotitesolidusdatarr.iloc[:,0].values,peridotitesolidusdatarr.iloc[:,1].values,linestyle=':',color='orange',label='Peridotite Solidus\n(Fiquet et al. 2010)')
    if 1:
        wcol = "#653700"
        wcol2 = "#FBDD7E"
        wcol3 = "#F97306"
        wop = 0.5
        wop2=0.75
        wop3=0.25
        ind = np.arange(0,11)
        print(stsmc.P[ind])
        axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],label='melt curve exp',alpha=wop,color=wcol)
        ind = np.arange(12,16)
        print(stsmc.P[ind])
        axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],alpha=wop,color=wcol)
        ind = np.arange(17,21)
        print(stsmc.P[ind])
        axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],alpha=wop,color=wcol)
        ind = np.arange(22,24)
        print(stsmc.P[ind])
        axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],alpha=wop2,color=wcol2,label='melt curve calc')
        ind = np.arange(25,28)
        print(stsmc.P[ind])
        axes[ai,aj].fill_between(stsmc.P[ind],stsmc.T_s[ind],stsmc.T_l[ind],alpha=wop2,color=wcol2)
        # Millot stv melt curve crossing
        axes[ai,aj].fill_between([470.,530.],[8000.,8000.],[8600.,8600.],alpha=0.5,color='blue',label='MM15')
 

    if 0:
        ind = np.where(stsmcsmooth.P <=Pt3)[0]
        axes[ai,aj].plot(stsmcsmooth.P[ind],stsmcsmooth.T[ind],'m-',linewidth=3,label='Proposed melt curve')
        ind = np.where(stsmcsmooth.P >=Pt3)[0]
        axes[ai,aj].plot(stsmcsmooth.P[ind],stsmcsmooth.T[ind],'m:',linewidth=3)#,label='Proposed melt curve')

    # density - pressure [3,0]
    ai=3
    aj=0
    
    #axes[ai,aj].scatter(qtz_hug_LA80.rho[6],qtz_hug_LA80.P[6],color='green',label=qtz_hug_LA80.label)
    #axes[ai,aj].plot(qtz_hug_KD.rho,qtz_hug_KD.P,color='green',label=qtz_hug_KD.label)

    # density - temperature [3,1]
    ai=3
    aj=1
    # Connolly JGR 2016 summary of critical points for silica
    axes[ai,aj].add_patch(Rectangle((0.5,5071), 0.5, 1000,  alpha=0.7, facecolor='lime', linewidth=0))
    axes[ai,aj].scatter(1900, -1.e-9, alpha=0.7, marker='s',color='lime', label='Connolly2016 CP')

    axes[ai,aj].plot(60.083/27.30,1999.,'s',color='red',label='Bourova Richet 1998')
    axes[ai,aj].plot(60.083/23.45,1700.,'d',color='orange',label='Bourova Richet 1998')
    axes[ai,aj].plot(60.083/23.51,1673.,'d',color='orange',label='')

    #axes[ai,aj].scatter(qtz_hug_LA80.rho[6],qtz_hug_LA80.T[6],color='green',label=qtz_hug_LA80.label)

    # density - pressure [3,0]
    ai=3
    aj=1    
    # Connolly JGR 2016
    axes[ai,aj].plot(C16bsrho,C16bstemp,'k-',label='C16 HSC-vWRK')
    axes[ai,aj].plot(C16bdrho,C16bdtemp,'k--',label='C16 HSC-CSRK')
    axes[ai,aj].plot(C16rsrho,C16rstemp,'r-',label='C16 SiO(l)-vWRK')
    # critical points
    axes[ai,aj].plot(1.036,6074,'ko')
    axes[ai,aj].plot(.555,5548,'ko')
    axes[ai,aj].plot(.782,5675,'ro')

axes[0,0].legend(fontsize=sfont-5)
axes[0,1].legend(fontsize=sfont-5)
axes[1,0].legend(fontsize=sfont-5)
axes[1,1].legend(fontsize=sfont-5)
axes[2,0].legend(fontsize=sfont-5)
axes[2,1].legend(fontsize=sfont-5)
axes[3,0].legend(fontsize=sfont-5)
axes[3,1].legend(fontsize=sfont-5)




plt.tight_layout(pad=0.6)
plt.savefig(MODELNAME+'-phase-diagram-fig.pdf', format='pdf', dpi=1000,transparent=True)
#--------

[13.7  14.5  20.   23.75 28.   34.   36.   37.5  40.   42.   46.  ]
[55. 60. 66. 74.]
[100. 102. 114. 117.]
[135. 155.]
[217.5 218.5   nan]

# 1-bar plots and 298 K plots

# Plot aneos vapor curves
sfont = 15
fig, axes = plt.subplots(2, 2, figsize=(14,14))
plt.subplots_adjust(wspace=0.4)

#--------
# temperature - entropy at 1 bar
ai=0
aj=0

axes[ai,aj].plot(NewEOS.onebar.T,NewEOS.onebar.S*1.E3,'+',color='blue',label='ANEOS 1 bar')

axes[ai,aj].plot(qtz_1bar_JANAF.T,qtz_1bar_JANAF.S,linestyle='-',color='brown',label=qtz_1bar_JANAF.label)
#axes[ai,aj].plot(janafgas_1bar.T,janafgas_1bar.S/1.e3,linestyle='-',color='orange',label='JANAF gas')
axes[ai,aj].scatter([3127,3127],[3.507,7.267],label='JANAF vap') 

axes[ai,aj].plot([3130,3130],[3.554,7.460],'kD',label='C16 1-bar boiling pt.')

axes[ai,aj].set_xlabel('Temperature (K)',size=sfont)
axes[ai,aj].set_ylabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_title('1-bar isobar',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(0,10)
axes[ai,aj].set_xlim(0,4000)


#----
# density-temperature at 1 bar
ai=0
aj=1

axes[ai,aj].plot(NewEOS.onebar.T,NewEOS.onebar.rho,'+',color='blue',label='ANEOS 1 bar')

axes[ai,aj].set_ylabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_xlabel('Temperature (K)',size=sfont)
axes[ai,aj].set_title('1-bar isobar',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(1,4)
axes[ai,aj].set_xlim(0,4000)

# plot data

axes[ai,aj].plot(1999.,60.083/27.30,'s',color='red',label='Bourova Richet 1998')
axes[ai,aj].plot(1700.,60.083/23.45,'d',color='orange',label='Bourova Richet 1998')
axes[ai,aj].plot(1673.,60.083/23.51,'d',color='orange',label='')

#----
# density-pressure at 298 K
ai=1
aj=0

iT0REF = np.where(NewEOS.T == T0REF)[0]
axes[ai,aj].plot(NewEOS.rho,NewEOS.P[iT0REF[0],:],'-',color='blue',label='ANEOS 298-K isotherm')
iT0REF = np.where(NewEOS.T == 2000)[0]
axes[ai,aj].plot(NewEOS.rho,NewEOS.P[iT0REF[0],:],'-.',color='red',label='ANEOS 2000-K isotherm')
iT0REF = np.where(NewEOS.T > 2600)[0]
axes[ai,aj].plot(NewEOS.rho,NewEOS.P[iT0REF[0],:],'-.',color='cyan',label='ANEOS 2600-K isotherm')
iT0REF = np.where(NewEOS.T == T0REF)[0]

# Finkelstein et al. 2014 Forsterite isotherm
K0T = 130.0 # +-9 GPa
K0Tp = 4.12 # +- 0.07
R0F = 3.22 # g/cm3
BM_pressure = 1.5*K0T*(np.power(NewEOS.rho/R0F,7./3.)-np.power(NewEOS.rho/R0F,5./3.)
                      )*(1.+0.75*(K0Tp-4.)*(np.power(NewEOS.rho/R0F,2./3.)-1.))
K0T = 130.0-9. # +-9 GPa
K0Tp = 4.12 # +- 0.07
BM_pressure_low = 1.5*K0T*(np.power(NewEOS.rho/R0F,7./3.)-np.power(NewEOS.rho/R0F,5./3.)
                      )*(1.+0.75*(K0Tp-4.)*(np.power(NewEOS.rho/R0F,2./3.)-1.))
K0T = 130.0+9. # +-9 GPa
K0Tp = 4.12 # +- 0.07
BM_pressure_high = 1.5*K0T*(np.power(NewEOS.rho/R0F,7./3.)-np.power(NewEOS.rho/R0F,5./3.)
                      )*(1.+0.75*(K0Tp-4.)*(np.power(NewEOS.rho/R0F,2./3.)-1.))
igood = np.where(BM_pressure < 48.)
axes[ai,aj].plot(NewEOS.rho[igood],BM_pressure[igood],'-',color='green',label='Fo 298-K isotherm\nFinkelstein et al. 2014')
igood = np.where(BM_pressure_low < 48.)
axes[ai,aj].plot(NewEOS.rho[igood],BM_pressure_low[igood],'--',color='green')
igood = np.where(BM_pressure_high < 48.)
axes[ai,aj].plot(NewEOS.rho[igood],BM_pressure_high[igood],'--',color='green')


# Jacobsen et al. 2008 MgO isotherm
K0T = 159.6 # +-0.6 GPa room temp isotherm
K0Tp = 3.74 # +-0.03 
R00 = (24.305+15.999)*4.*1.6726/(74.697*1.E-3)/1.E3 # kg/m3->g/cm3 
print(R00)
BM_pressure = 1.5*K0T*(np.power(NewEOS.rho/R00,7./3.)-np.power(NewEOS.rho/R00,5./3.)
                      )*(1.+0.75*(K0Tp-4.)*(np.power(NewEOS.rho/R00,2./3.)-1.))
igood = np.where((BM_pressure > 25.)&(BM_pressure < 125.))
axes[ai,aj].plot(NewEOS.rho[igood],BM_pressure[igood],'-',color='brown',label='MgO 298-K isotherm\nJacobsen et al. 2008')

# Vanpeteghem et al. 2006 silicate perovskite isotherm
K0T = 253. # +-1 GPa room temp isotherm
K0Tp = 4.0 
R00 = (24.305+28.085+3.*15.999)*4.*1.6726/(162.51*1.E-3)/1.E3 # kg/m3->g/cm3 
print(R00)
BM_pressure = 1.5*K0T*(np.power(NewEOS.rho/R00,7./3.)-np.power(NewEOS.rho/R00,5./3.)
                      )*(1.+0.75*(K0Tp-4.)*(np.power(NewEOS.rho/R00,2./3.)-1.))
igood = np.where((BM_pressure > 25.)&(BM_pressure < 125.))
axes[ai,aj].plot(NewEOS.rho[igood],BM_pressure[igood],'-',color='orange',label='MgSiO$_3$ Pv 298-K isotherm\nVanpeteghem et al. 2006')

# Hosono et al. nat Geo hard sphere model isotherm 2000 K
K0T = 25. # GPa
K0Tp = 5.2
R00 = 2.65
BM_pressure = 1.5*K0T*(np.power(NewEOS.rho/R00,7./3.)-np.power(NewEOS.rho/R00,5./3.)
                      )*(1.+0.75*(K0Tp-4.)*(np.power(NewEOS.rho/R00,2./3.)-1.))
igood = np.where(BM_pressure < 40.)
axes[ai,aj].plot(NewEOS.rho[igood],BM_pressure[igood],'-',color='red',label='Karato 2000K hard sphere isotherm')

axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_title('298-K isotherm',size=sfont)
axes[ai,aj].set_xlim(2.5,6)
axes[ai,aj].set_ylim(0,140)

#----
if 0:
    # density-pressure at 298 K
    ai=1
    aj=1

    axes[ai,aj].plot(mantle.density,mantle.pressure,'-',color='blue',label='ANEOS '+str(np.round(np.interp(1.E-4,mantle.pressure,mantle.temperature)))+' K isentrope',markersize=10)
    axes[ai,aj].plot(PREM.density/1.E3,PREM.pressure,'-',color='orange',label='PREM')

    mantle = isentrope_class() # code in eostable.py
    Tisentrope = 1000. #K
    Smantle = np.interp(Tisentrope,NewEOS.onebar.T,NewEOS.onebar.S)

    # loop across all densities and extract the values for the requested isentrope
    for i in range(0,NewEOS.ND):
        ind = np.where((NewEOS.S[:,i] > 0))[0]
        interpfunction = interpolate.interp1d(NewEOS.S[ind,i],NewEOS.P[ind,i]) # MJ/K/kg, GPa
        mantle.pressure = np.append(mantle.pressure,interpfunction(Smantle)) # GPa
        interpfunction = interpolate.interp1d(NewEOS.S[ind,i],NewEOS.T[ind]) # MJ/K/kg, GPa
        mantle.temperature = np.append(mantle.temperature,interpfunction(Smantle)) # GPa
    mantle.density = NewEOS.rho # g/cm3

    print('Calculated isentrope with ',np.round(np.interp(1.E-4,mantle.pressure,mantle.temperature)),' K potential temperature')
    axes[ai,aj].plot(mantle.density,mantle.pressure,'-.',color='blue',label='ANEOS '+str(np.round(np.interp(1.E-4,mantle.pressure,mantle.temperature)))+' K isentrope',markersize=10)

    axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
    axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
    axes[ai,aj].tick_params(labelsize=sfont)
    axes[ai,aj].set_title('PREM and '+str(np.round(np.interp(1.E-4,mantle.pressure,mantle.temperature)))+' K ANEOS model isentrope',size=sfont)
    axes[ai,aj].set_xlim(2,6)
    axes[ai,aj].set_ylim(0,140)

# don't show a plot in lower right
#axes[1,1].axis("off")


axes[0,0].legend(fontsize=sfont-5)
axes[0,1].legend(fontsize=sfont-5)
axes[1,0].legend(fontsize=sfont-5)
axes[1,1].legend(fontsize=sfont-5)

3.6099158145574792
4.132848346563289

/var/folders/4z/xws_ywgj78sc494vpwymn9bw0000gn/T/ipykernel_87866/2785270329.py:159: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
  axes[1,1].legend(fontsize=sfont-5)

<matplotlib.legend.Legend at 0x163cda710>

# plot the gamma and isotherms

#NewEOS.R0REF=4.0
#NewEOS.gamma0=1.1

# ANEOS Formulation for Gamma and theta
theta0=abs(NewEOS.theta0)
rhoarr = (1+np.arange(1000))/1000.*NewEOS.R0REF*10.
gammaarr = (NewEOS.gamma0*NewEOS.R0REF/rhoarr
           )+NewEOS.C24*np.power(1.0-NewEOS.R0REF/rhoarr,2.)*(1.-NewEOS.C60)+(NewEOS.gamma0 + (NewEOS.C24-NewEOS.gamma0)*np.power(1.-NewEOS.R0REF/rhoarr,2.))*NewEOS.C60
s3 = NewEOS.C61+1.
C14 = theta0 * np.exp(1.5*s3-(2.-0.5*NewEOS.C60)*NewEOS.gamma0)/(np.power(NewEOS.R0REF,s3))
C16 = (s3-(2.-NewEOS.C60)*NewEOS.gamma0)/(NewEOS.R0REF*NewEOS.R0REF)
C17 = ((3.-NewEOS.C60)*NewEOS.gamma0-2.*s3)/NewEOS.R0REF
gammaarrlow = C16*rhoarr*rhoarr+C17*rhoarr+s3
thetaarr = theta0*np.exp((1.-NewEOS.C60)*(1.-NewEOS.R0REF/rhoarr)*NewEOS.gamma0-0.5*(NewEOS.R0REF/rhoarr)*(3.-(NewEOS.R0REF/rhoarr)*(4-(NewEOS.R0REF/rhoarr))))*np.power((rhoarr/NewEOS.R0REF),NewEOS.C24)
thetalowarr = C14*np.exp(rhoarr*(C17+0.5*C16*rhoarr))*np.power(rhoarr,s3)


#-------- plots
sfont = 15
fig, axes = plt.subplots(1, 2, figsize=(14,8))
plt.subplots_adjust(wspace=0.4)


# -----
ai=0
indlow = np.where(rhoarr <= NewEOS.R0REF)
indhigh = np.where(rhoarr >= NewEOS.R0REF)
axes[ai].plot(rhoarr[indhigh[0]],gammaarr[indhigh[0]],label="ANEOS Gamma-high",color='black')
axes[ai].plot(rhoarr[indlow[0]], gammaarrlow[indlow[0]],label="ANEOS Gamma-low",color='black',ls='--')

axes[ai].plot(NewEOS.R0REF,NewEOS.gamma0,'+',label="ANEOS Gamma0",color='black')

axes[ai].set_ylim(0,2.5)
axes[ai].set_xlim(0,8)

axes[ai].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai].set_ylabel('Gamma (-)',size=sfont)
axes[ai].tick_params(labelsize=sfont)

# solid forsterite at STP
#axes[ai].plot([3.22],[1.29],'s',color='red',label='Fo solid STP')

# solid alpha-quartz at STP Rammo & Farid Powder Diffraction 1994 
axes[ai].plot([2.65],[0.641],'s',color='red',label='a-qtz solid STP')

# Thomas and Asimow 2013 liquid forsterite shock data
gamma_0=0.396
q=-2.02
q_u=1.03
Rho_0=2597. # kg/m^3
asi_rhot = np.array([3739.79122056, 3335.90237636,4691.02371542 ,4088.40354767, 4347.15217391 ,4681.1865942 ])
##asi_rhot=np.array([3340.0, 3740.0, 4100.0, 4680.0, 5290.0, 5420.0]) # BAD TABULATED VALUES in TA13
asi_gammat=gamma_0*np.power(Rho_0/asi_rhot,q)

asi_rho = np.append([2597.],asi_rhot)/1.e3
asi_gamma = np.append([0.396],asi_gammat)

#print(asi_rho)
#print(asi_gamma)
#[3739.79122056 3335.90237636 4691.02371542 4088.40354767 4347.15217391 4681.1865942 ]
#[0.82720568 0.65667847 1.30744208 0.99037675 1.12107749 1.30190973]

# Erik's error bars monte carlo'd for the asimow data
#gamma [ 0.82271901  0.65485213  1.29584059  0.97276836  1.11171104  1.29195917]
asi_gamma_sig_low=np.asarray([0., 0.25715384 , 0.14948592 , 0.5898531,   0.36154228 , 0.45607297 , 0.58790934])
asi_gamma_sig_high=np.asarray([0., 0.37855329,  0.19525818 , 1.09577441,  0.6137275 ,  0.78965792 , 1.08676142])

#axes[ai].plot(asi_rho,asi_gamma,'o',color='orange',label='Forsterite liq.\nTA13 corrected')
#for i in range(0,len(asi_gamma)):
#    axes[ai].plot([asi_rho[i],asi_rho[i]],[asi_gamma[i]-asi_gamma_sig_low[i],asi_gamma[i]+asi_gamma_sig_high[i]],linestyle='-',color='orange')

# these values from the submitted paper
#axes[ai].errorbar(SR_rho,SR_gamma,yerr=SR_gamma_e, xerr=SR_rho_e, fmt='g+', label ='Forsterite liq. Z data\nDavies et al. 2019')

gamma_0=0.5
Rho_0=2.65
Karato_gammar=np.arange(.7,1,.01)
q=-2
Karato_gammat=gamma_0*np.power(Karato_gammar,q)
#axes[ai].plot(Rho_0/(Karato_gammar),Karato_gammat,label='Karato gamma model')
#print(Rho_0/Karato_gammar)
#print(Karato_gammat)

axes[ai].plot(qtz_hug_McCoy16.rho,qtz_hug_McCoy16.gamma,'+',color='red',label=qtz_hug_McCoy16.label)
axes[ai].plot(karkigcvorig.rho,karkigcvorig.gamma,'D',label='Karki 2007')
# -----
ai=1
axes[ai].plot(rhoarr[indhigh[0]],thetaarr[indhigh[0]],label="ANEOS Theta-high",color='black')
axes[ai].plot(rhoarr[indlow[0]],thetalowarr[indlow[0]],label="ANEOS Theta-low",color='black',ls='--')
axes[ai].plot(NewEOS.R0REF,theta0,'+',label="Theta0",color='black')

axes[ai].set_ylim(0,10000)
axes[ai].set_xlim(0,8)

axes[ai].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai].set_ylabel('Theta (K)',size=sfont)
axes[ai].tick_params(labelsize=sfont)


axes[0].legend(fontsize=sfont-5)
axes[1].legend(fontsize=sfont-5)

<matplotlib.legend.Legend at 0x318db5e50>

sfont = 15
fig, axes = plt.subplots(2, 2, figsize=(14,14))
plt.subplots_adjust(wspace=0.4)


#-----------------------------
# Pick isotherms and set color scale
ixtemparr = [298.,1000.,2050.,2300.,5000.,8000.,20000.]
# set color scale same for all plots
tcolmin=0.
tcolmax=20000.

# find the index closest to the selected temperatures -- no interpolation
ixindarr=np.zeros(len(ixtemparr))
ixlabelarr = []
for i in np.arange(len(ixtemparr)):
    tmpi=int(np.round(np.interp(ixtemparr[i],NewEOS.T,np.arange(NewEOS.NT))))
    #print(tmpi,antarr[tmpi])
    ixindarr[i] = tmpi
    ixlabelarr.append(str(ixtemparr[i])+' K')


#-----------------------
ai=0
aj=0
for i in np.arange(len(ixtemparr)):
    icolscale = (ixtemparr[i]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.cv[int(ixindarr[i]),:]*1000.,c=col,linestyle='-',label=ixlabelarr[i])

#axes[ai,aj].plot(QMDH_r1,QMDH_cv*1000.,styleqmd[0],markersize=styleqmd[2],color=styleqmd[1],label='QMD Hugoniot\n(Root et al. 2018')
#axes[ai,aj].plot(2.597,1.73736,'^',color='orange',label='1 bar 2273 K liquid Cv\n(Thomas & Asimow 20133)')
#axes[ai,aj].plot(3.,180./(NewEOS.FMW/1.E3)/1.E3,'s',color='orange',label='1 bar 2000 K solid Cv\n(Gillet et al. 1991)')
axes[ai,aj].set_ylim(1.,5.)
axes[ai,aj].set_xlim(0,10)

axes[ai,aj].set_ylabel('Specific heat capacity Cv (kJ/K/kg)',size=sfont)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-----------------------
ai=0
aj=1
for i in np.arange(len(ixtemparr)):
    icolscale = (ixtemparr[i]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.cv[int(ixindarr[i]),:]*1000.*1000./7/8.314*0.14,c=col,linestyle='-',label=ixlabelarr[i])


#axes[ai,aj].plot(QMDH_r1,QMDH_cv*1.E6/7/8.314*0.14,styleqmd[0],markersize=styleqmd[2],color=styleqmd[1],label='QMD Hugoniot\n(Root et al. 2018')
#axes[ai,aj].plot(2.597,1737./7/8.314*0.140,'^',color='orange',label='1 bar 2273 K liquid Cv\n(Thomas & Asimow 2013)')
#axes[ai,aj].plot(3.,180./(NewEOS.FMW/1.E3)/7/8.314*0.140,'s',color='orange',label='1 bar 2000 K solid Cv\n(Gillet et al. 1991)')
axes[ai,aj].set_ylim(1.,10.)
axes[ai,aj].set_xlim(0,10)

axes[ai,aj].set_ylabel('Specific heat capacity Cv (nR)',size=sfont)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-----------------------
# plot sound speed along isotherms
ai=1
aj=0
for i in np.arange(len(ixtemparr)):
    icolscale = (ixtemparr[i]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.P[int(ixindarr[i]),:],NewEOS.cs[int(ixindarr[i]),:]/1.e5,c=col,linestyle='-',label=ixlabelarr[i])

#axes[ai,aj].plot(1.e-4,2.674,'d',color='orange',label='Thomas & Asimow 2013 2273 K')

axes[ai,aj].set_ylim(0,10)
axes[ai,aj].set_xlim(-10,40)

axes[ai,aj].set_ylabel('Bulk Sound Speed (km/s)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


#-----------------------
# plot sound speed along isotherms
ai=1
aj=1
for i in np.arange(len(ixtemparr)):
    icolscale = (ixtemparr[i]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.P[int(ixindarr[i]),:],NewEOS.cs[int(ixindarr[i]),:]/1.e5,c=col,linestyle='-',label=ixlabelarr[i])

#axes[ai,aj].plot(1.e-4,2.674,'d',color='orange',label='Thomas & Asimow 2013 2273 K')

axes[ai,aj].plot(qtz_hug_McCoy16.P,qtz_hug_McCoy16.cs,'+',color='red',label=qtz_hug_McCoy16.label)


axes[ai,aj].set_ylim(0,20)
axes[ai,aj].set_xlim(0,700)

axes[ai,aj].set_ylabel('Bulk Sound Speed (km/s)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

axes[0,0].legend(fontsize=sfont-5)
axes[0,1].legend(fontsize=sfont-5)
axes[1,0].legend(fontsize=sfont-5)
axes[1,1].legend(fontsize=sfont-5)

<matplotlib.legend.Legend at 0x318dc0dd0>

# Load older ANEOS models for comparison
# Load in Melosh 2013 model
# READ IN NEW ANEOS MODEL and fill the extEOStable class object
# source in eostable.py
#------------------------------------------------------------------
ORIGMODELNAME = MODELNAME
subdir='silica-Melosh2007/'
MODELNAME = 'Silica_Melosh_2007'
# Header information must all be compatible with float format
MATID = 1.0        # MATID number
DATE = 070513.     # Date as a single 6-digit number YYMMDD
VERSION = 1.0      # ANEOS Parameters Version number
FMN = 30.          # Formula weight in atomic numbers for silica
FMW = 60.083       # Formula molecular weight (g/mol)
# The following define the default initial state for material in the 201 table
R0REF   = 2.65     # g/cm3 *** R0REF is inserted into the density array
K0REF   = 3.58E11  # dynes/cm2
T0REF   = 298.     # K -- *** T0REF is inserted into the temperature array
P0REF   = 1.E6     # dynes/cm2 -- this defines the principal Hugoniot calculated below

MeloshEOS  = etab.extEOStable() # FIRST make new empty EOS object
MeloshEOS.loadextsesame(subdir+'NEW-SESAME-EXT.TXT') # LOAD THE EXTENDED 301 SESAME FILE GENERATED BY STSM VERSION OF ANEOS
MeloshEOS.loadstdsesame(subdir+'NEW-SESAME-STD.TXT') # LOAD THE STANDARD 301 SESAME FILE GENERATED BY STSM VERSION OF ANEOS
MeloshEOS.MODELNAME = MODELNAME # string set above in user input
MeloshEOS.MDQ = np.zeros((MeloshEOS.NT,MeloshEOS.ND)) # makes the empty MDQ array
#print(MeloshEOS.units) # these are the default units for SESAME rho-T tables
#'Units: g/cm3, K, GPa, MJ/kg, MJ/kg, MJ/K/kg, cm/s, MJ/K/kg, KPA flag. 2D arrays are (NT,ND).'

# Add the header info to the table. This could be done during the loading. 
# if made from this notebook, these values are set in the user-input above.
# ** MAKE SURE THEY MATCH ANEOS.INPUT **
MeloshEOS.MATID   = MATID
MeloshEOS.DATE    = DATE
MeloshEOS.VERSION = VERSION
MeloshEOS.FMN     = FMN
MeloshEOS.FMW     = FMW
MeloshEOS.R0REF   = R0REF
MeloshEOS.K0REF   = K0REF
MeloshEOS.T0REF   = T0REF
MeloshEOS.P0REF   = P0REF
#
# Load the information from ANEOS.INPUT and ANEOS.OUTPUT
MeloshEOS.loadaneos(aneosinfname=subdir+'ANEOS.INPUT',aneosoutfname=subdir+'ANEOS.OUTPUT',silent=True)
#
# calculate the principal Hugoniot from the table
MeloshEOS.calchugoniot(r0=MeloshEOS.R0REF,t0=MeloshEOS.T0REF)
# Calculate the fused silica Hugoniot from the table
MeloshEOS.calcOffEOSHugoniot(r0=2.20,t0=298.,p0=1.e-4,e0=0.,r1=2.65)
#
# calculate the 1-bar profile; loop over temp
MeloshEOS.onebar.T = np.zeros(MeloshEOS.NT)
MeloshEOS.onebar.S = np.zeros(MeloshEOS.NT)
MeloshEOS.onebar.rho = np.zeros(MeloshEOS.NT)
it0 = np.where(MeloshEOS.T >= MeloshEOS.T0REF)[0]
id0 = np.arange(MeloshEOS.ND)#np.where(MeloshEOS.rho >= 0.8*MeloshEOS.R0REF)[0]
for iit in range(0,MeloshEOS.NT):
    MeloshEOS.onebar.T[iit] = MeloshEOS.T[iit]
    MeloshEOS.onebar.S[iit] = np.interp(1.E-4,MeloshEOS.P[iit,id0],MeloshEOS.S[iit,id0])
    MeloshEOS.onebar.rho[iit] = np.interp(1.E-4,MeloshEOS.P[iit,id0],MeloshEOS.rho[id0])
    #print(iit,MeloshEOS.onebar.rho[iit],MeloshEOS.onebar.T[iit],MeloshEOS.onebar.S[iit]*1.E3)

rho0= 2.65
r1= 2.65
367 22 2.2 298.0 0.0001 0.0 0.0 3.67465878

# Display the initial conditions for the Hugoniot in ANEOS.OUTPUT    
print('\nRequested ANEOS Hugoniot initial state from (THUG,RHUG) (red curves in plots below):')
print('Initial Pressure (GPa):',NewEOS.anhug.P[0])
print('Initial Temperature (K):',NewEOS.anhug.T[0])
print('Initial Sp. Entropy (MJ/K/kg):',NewEOS.anhug.S[0])
print('Initial Sp. Internal Energy (MJ/kg):',NewEOS.anhug.U[0])
print('Initial Density (g/cm3):',NewEOS.anhug.rho[0])
#

#========== PLOTTING STYLES FOR LOTS OF PLOTS =================
styleqmd = ['^','red',5,'Fo QMD calc (Root et al. 2018)']
stylezdat = ['^','black',5,'Fo Z data (Root et al. 2018)']
stylegg = ['+','green',5,'Fo Gas gun data\n(Mosenfelder et al. 2007)']
styleliq = ['^','orange',7,'Liquid Fo data corrected\n(Thomas & Asimow 2013)']
styleliqhug = ['-','orange',1,'SLVT ANEOS Liquid 2273K Hugoniot']
stylezfit = ['-','green',1,'Fo Z fit (Davies et al. 2019)']
styleanhug = ['-','red',1,'SLVT ANEOS Glass-tension Hugoniot']
stylestphug = ['-','green',3,'SLVT ANEOS Qtz Hugoniot']
styleanphase = ['-','black',1,'SLVT ANEOS Phase Boundaries']
stylephasept = ['o','black',6]

styleanFS = ['--','blue',3,'SLVT ANEOS FS Hugoniot']

styleqtzhugmelosh = ['-','goldenrod',2,'M07 ANEOS Qtz Hugoniot']
stylefshugmelosh = [':','goldenrod',2,'M07 ANEOS FS Hugoniot']
styleanphasemelosh = ['-','goldenrod',1,'M07 ANEOS Phase Boundaries']


#======================LOTS OF PLOTS============================
# plot the Hugoniots for comparison
sfont = 15
fig, axes = plt.subplots(nrows=4, ncols=2, figsize=(14,22))
plt.subplots_adjust(wspace=0.4)

#-------------
# Pressure-Density Data vs. ANEOS
ai=0
aj=0
#axes[ai,aj].plot(NewEOS.anhug.rho,NewEOS.anhug.P,styleanhug[0],linewidth=styleanhug[2],color=styleanhug[1],label=styleanhug[3])

axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.rho,NewEOS.hugo.P,styleanFS[0],linewidth=styleanFS[2],color=styleanFS[1],label=styleanFS[3])

axes[ai,aj].plot(MeloshEOS.hug.rho,MeloshEOS.hug.P,styleqtzhugmelosh[0],linewidth=styleqtzhugmelosh[2],color=styleqtzhugmelosh[1],label=styleqtzhugmelosh[3])
axes[ai,aj].plot(MeloshEOS.hugo.rho,MeloshEOS.hugo.P,stylefshugmelosh[0],linewidth=stylefshugmelosh[2],color=stylefshugmelosh[1],label=stylefshugmelosh[3])
     
if 0:
    axes[ai,aj].scatter(qtz_hug_W62.rho,qtz_hug_W62.P,marker='+',color='purple',label=qtz_hug_W62.label)
    axes[ai,aj].scatter(qtz_hug_L3.rho,qtz_hug_L3.P,marker='+',color='blue',label=qtz_hug_L3.label)
    axes[ai,aj].plot(qtz_hug_KD.rho,qtz_hug_KD.P,color='green',linewidth=5,label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.rho,fs_hug_RTK19.P,':',color='green',linewidth=5,label=fs_hug_RTK19.label)

if 1:
    axes[ai,aj].errorbar(shockdat_list[K13qtz['id']].rho,shockdat_list[K13qtz['id']].P,fmt='gd',markersize=3,\
             xerr=shockdat_list[K13qtz['id']].rho_std,yerr=shockdat_list[K13qtz['id']].P_std,label=K13qtz['lab'])
    ind = np.where(fs_hug_RTK19.P > 74.)[0]
    axes[ai,aj].plot(fs_hug_RTK19.rho[ind],fs_hug_RTK19.P[ind],'b--',label='R19 FS',linewidth=1)

axes[ai,aj].set_ylim(0,1200)
axes[ai,aj].set_xlim(4.,8)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Pressure - Temperature Data vs. ANEOS
ai=0
aj=1

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.T/1.e3,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.P,NewEOS.hugo.T/1.e3,styleanFS[0],linewidth=styleanFS[2],color=styleanFS[1],label=styleanFS[3])

axes[ai,aj].plot(MeloshEOS.hug.P,MeloshEOS.hug.T/1.e3,styleqtzhugmelosh[0],linewidth=styleqtzhugmelosh[2],color=styleqtzhugmelosh[1],label=styleqtzhugmelosh[3])
axes[ai,aj].plot(MeloshEOS.hugo.P,MeloshEOS.hugo.T/1.e3,stylefshugmelosh[0],linewidth=stylefshugmelosh[2],color=stylefshugmelosh[1],label=stylefshugmelosh[3])

if 0:
    axes[ai,aj].plot(qtz_hug_KD.P,qtz_hug_KD.T/1.e3,color='green',linewidth=5,label='Qtz Brygoo et al. 2015')
    axes[ai,aj].plot(fs_hug_Hicks06.P,fs_hug_Hicks06.T/1.e3,':',color='green',linewidth=5,label=fs_hug_Hicks06.label)

    axes[ai,aj].scatter(qtz_hug_LA80.P,qtz_hug_LA80.T/1.E3,marker='+',color='green',label='Qtz Lyzenga & Ahrens 1980')
    axes[ai,aj].plot(stv_mc_MM15.Pl,stv_mc_MM15.T/1.E3,color='black',label=stv_mc_MM15.label)
if 1:
    axes[ai,aj].errorbar(fs_hug_Hicks06.P,fs_hug_Hicks06.T/1.e3,fmt='b+:',\
             yerr=fs_hug_Hicks06.T_std/1.e3,label='H06 FS')
    axes[ai,aj].plot(qtz_hug_Hicks06.P,qtz_hug_Hicks06.T/1.e3,'g+:',label='H06 qtz')

axes[ai,aj].set_xlim(0,1200)
axes[ai,aj].set_ylim(0,70)
axes[ai,aj].set_ylabel('Temperature (1000 K)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Particle Velocity vs. Shock Velocity; Data vs. ANEOS
ai=1
aj=0

axes[ai,aj].plot(NewEOS.hug.up,NewEOS.hug.us,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.up,NewEOS.hugo.us,styleanFS[0],linewidth=styleanFS[2],color=styleanFS[1],label=styleanFS[3])

axes[ai,aj].plot(MeloshEOS.hug.up,MeloshEOS.hug.us,styleqtzhugmelosh[0],linewidth=styleqtzhugmelosh[2],color=styleqtzhugmelosh[1],label=styleqtzhugmelosh[3])
axes[ai,aj].plot(MeloshEOS.hugo.up,MeloshEOS.hugo.us,stylefshugmelosh[0],linewidth=stylefshugmelosh[2],color=stylefshugmelosh[1],label=stylefshugmelosh[3])

if 0:
    axes[ai,aj].scatter(qtz_hug_W62.up,qtz_hug_W62.us,marker='+',color='purple',label=qtz_hug_W62.label)
    axes[ai,aj].scatter(qtz_hug_L3.up,qtz_hug_L3.us,marker='+',color='blue',label=qtz_hug_L3.label)
    axes[ai,aj].plot(qtz_hug_KD.up,qtz_hug_KD.us,color='green',label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.up,fs_hug_RTK19.us,':',color='green',label=fs_hug_RTK19.label)

if 1:
    axes[ai,aj].errorbar(shockdat_list[K13qtz['id']].up,shockdat_list[K13qtz['id']].us,fmt='gd',markersize=3,\
             xerr=shockdat_list[K13qtz['id']].up_std,yerr=shockdat_list[K13qtz['id']].us_std,label=K13qtz['lab'])
    ind = np.where(fs_hug_RTK19.P > 74.)[0]
    axes[ai,aj].plot(fs_hug_RTK19.up[ind],fs_hug_RTK19.us[ind],'b--',label='R19 FS',linewidth=1)


axes[ai,aj].set_xlim(0,15)
axes[ai,aj].set_ylim(0,25)
axes[ai,aj].set_xlabel('Particle Velocity (km/s)',size=sfont)
axes[ai,aj].set_ylabel('Shock Velocity (km/s)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Sp. Entropy vs Pressure; Calculated from experiments vs. ANEOS
ai=1
aj=1

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.S*1.e3,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.P,NewEOS.hugo.S*1.e3,styleanFS[0],linewidth=styleanFS[2],color=styleanFS[1],label=styleanFS[3])

axes[ai,aj].plot(MeloshEOS.hug.P,MeloshEOS.hug.S*1.e3,styleqtzhugmelosh[0],linewidth=styleqtzhugmelosh[2],color=styleqtzhugmelosh[1],label=styleqtzhugmelosh[3])
axes[ai,aj].plot(MeloshEOS.hugo.P,MeloshEOS.hugo.S*1.e3,stylefshugmelosh[0],linewidth=stylefshugmelosh[2],color=stylefshugmelosh[1],label=stylefshugmelosh[3])

if 0:
    # Kraus entropy fit only good to 800 GPa
    ind = np.where(qtz_hug_KD.P <= 800.)
    axes[ai,aj].plot(qtz_hug_KD.P[ind],qtz_hug_KD.S[ind],color='green',label='Qtz Kraus et al. 2012')
    Kraus_serr = (127.3+0.328*qtz_hug_KD.P[ind]-0.000193*np.power(qtz_hug_KD.P[ind],2))/1.e3 # kJ/K/kg
    axes[ai,aj].fill_between(qtz_hug_KD.P[ind],(qtz_hug_KD.S[ind]-Kraus_serr),(qtz_hug_KD.S[ind]+Kraus_serr),color=stylezfit[1],alpha=0.2)

# silicadb
if 1:
    entropyshift = 220.
    penter = 111.
    kqtzparr = np.arange(100,500) # GPa
    kqtzsarr = 1793+22.4*kqtzparr-0.0426*np.power(kqtzparr,2)+4.319e-05*np.power(kqtzparr,3)-1.742e-8*np.power(kqtzparr,4)
    kqtzserrarr = 127.3+0.328*kqtzparr-0.000193*np.power(kqtzparr,2)
    ind = np.where(kqtzparr >= penter)[0]
    axes[ai,aj].plot(kqtzparr[ind],kqtzsarr[ind]/1.e3,'g-',linewidth=1,label='K12 qtz Hugoniot-liq region')
    axes[ai,aj].plot(kqtzparr[ind],(kqtzsarr[ind]-entropyshift)/1.e3,'g-',linewidth=3,label='K12 Qtz Hugoniot shift')
    axes[ai,aj].fill_between(kqtzparr[ind],(kqtzsarr[ind]-entropyshift-kqtzserrarr[ind])/1.e3,(kqtzsarr[ind]-entropyshift+kqtzserrarr[ind])/1.e3,color=stylezfit[1],alpha=0.2)
   
axes[ai,aj].set_ylim(0,8.)
axes[ai,aj].set_xlim(0,1000)
axes[ai,aj].set_ylabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Pressure-Density Data vs. ANEOS   LOW PRESSURE + LIQUID HUGONIOT
ai=2
aj=0
# ANEOS REFERENCE POINT (STP) Hugoniots
axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.rho,NewEOS.hugo.P,styleanFS[0],linewidth=styleanFS[2],color=styleanFS[1],label=styleanFS[3])

axes[ai,aj].plot(MeloshEOS.hug.rho,MeloshEOS.hug.P,styleqtzhugmelosh[0],linewidth=styleqtzhugmelosh[2],color=styleqtzhugmelosh[1],label=styleqtzhugmelosh[3])
axes[ai,aj].plot(MeloshEOS.hugo.rho,MeloshEOS.hugo.P,stylefshugmelosh[0],linewidth=stylefshugmelosh[2],color=stylefshugmelosh[1],label=stylefshugmelosh[3])

if 0:
    axes[ai,aj].scatter(qtz_hug_W62.rho,qtz_hug_W62.P,marker='+',color='purple',label=qtz_hug_W62.label)
    axes[ai,aj].scatter(qtz_hug_L3.rho,qtz_hug_L3.P,marker='+',color='blue',label=qtz_hug_L3.label)
    axes[ai,aj].plot(qtz_hug_KD.rho,qtz_hug_KD.P,color='green',label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.rho,fs_hug_RTK19.P,':',color='green',label=fs_hug_RTK19.label)

# silicadb
if 1:
    axes[ai,aj].errorbar(shockdat_list[L83qtzhug['id']].rho,shockdat_list[L83qtzhug['id']].P,fmt='go',mfc='white',\
                 yerr=shockdat_list[L83qtzhug['id']].P_std,label=L83qtzhug['lab'])
    axes[ai,aj].errorbar(shockdat_list[W62qtz['id']].rho,shockdat_list[W62qtz['id']].P,fmt='gx',label=W62qtz['lab'])
    #plt.errorbar(shockdat_list[A23qtz['id']].rho,shockdat_list[A23qtz['id']].P,fmt='gD',mfc='white',label=A23qtz['lab'],\
    #            xerr=shockdat_list[A23qtz['id']].rho_std,yerr=shockdat_list[A23qtz['id']].P_std)
    axes[ai,aj].errorbar(shockdat_list[T20qtz['id']].rho,shockdat_list[T20qtz['id']].P,fmt='<',mfc='white',color='grey',\
                 xerr = shockdat_list[T20qtz['id']].rho_std, yerr=shockdat_list[T20qtz['id']].P_std, label=T20qtz['lab'])
    #print(shockdat_list[T20qtz['id']].up,shockdat_list[T20qtz['id']].us)
    #res = np.polyfit(shockdat_list[T20qtz['id']].up,shockdat_list[T20qtz['id']].us,1)
    #print(res)
    axes[ai,aj].errorbar(shockdat_list[L04qtz['id']].rho,shockdat_list[L04qtz['id']].P,fmt='gv',mfc='white',\
             xerr=shockdat_list[L04qtz['id']].rho_std,yerr=shockdat_list[L04qtz['id']].P_std,label=L04qtz['lab'])

    ind = np.where(qtz_hug_KD.P>112)[0]
    axes[ai,aj].plot(qtz_hug_KD.rho[ind],qtz_hug_KD.P[ind],'g-',linewidth=3,label='KD90 liq fit')
    axes[ai,aj].errorbar(qtzdf['rho'][indqtz],qtzdf['P'][indqtz],fmt='+',color='grey',label='M80 qtz')

axes[ai,aj].set_ylim(-5,200)
axes[ai,aj].set_xlim(2.5,6.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)

#-------------
# Pressure-Density Data vs. ANEOS   EXTRAPOLATED HIGH PRESSURE
ai=2
aj=1
axes[ai,aj].plot(NewEOS.hug.rho,NewEOS.hug.P,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hugo.rho,NewEOS.hugo.P,styleanFS[0],linewidth=styleanFS[2],color=styleanFS[1],label=styleanFS[3])

axes[ai,aj].plot(MeloshEOS.hug.rho,MeloshEOS.hug.P,styleqtzhugmelosh[0],linewidth=styleqtzhugmelosh[2],color=styleqtzhugmelosh[1],label=styleqtzhugmelosh[3])
axes[ai,aj].plot(MeloshEOS.hugo.rho,MeloshEOS.hugo.P,stylefshugmelosh[0],linewidth=stylefshugmelosh[2],color=stylefshugmelosh[1],label=stylefshugmelosh[3])

if 0:
    axes[ai,aj].plot(qtz_hug_KD.rho,qtz_hug_KD.P,color='green',label=qtz_hug_KD.label)
    axes[ai,aj].plot(fs_hug_RTK19.rho,fs_hug_RTK19.P,':',color='green',label=fs_hug_RTK19.label)

# silicadb
if 1:
    axes[ai,aj].errorbar(shockdat_list[W62fs['id']].rho,shockdat_list[W62fs['id']].P,fmt='bx',label=W62fs['lab'])
    axes[ai,aj].errorbar(shockdat_list[L83fshug['id']].rho,shockdat_list[L83fshug['id']].P,fmt='bo',label=L83fshug['lab'])
    axes[ai,aj].errorbar(shockdat_list[B19fs['id']].rho,shockdat_list[B19fs['id']].P,fmt='b>',xerr=shockdat_list[B19fs['id']].rho_std,\
                 yerr=shockdat_list[B19fs['id']].P_std,mfc='white',label=B19fs['lab'])
    axes[ai,aj].errorbar(shockdat_list[T18fs['id']].rho,shockdat_list[T18fs['id']].P,fmt='b<',mfc='white',\
                 xerr=shockdat_list[T18fs['id']].rho_std,yerr=shockdat_list[T18fs['id']].P_std,label=T18fs['lab'])

    #plt.errorbar(np.asarray(shockdat_list[A23fs['id']].rho)/1000.,shockdat_list[A23fs['id']].P,fmt='bD',mfc='white',label=A23fs['lab'])
    #plt.errorbar(np.asarray(shockdat_list[A23fs['id']].rho),shockdat_list[A23fs['id']].P,fmt='bD',mfc='white',label=A23fs['lab'],\
    #            xerr=shockdat_list[A23fs['id']].rho_std,yerr=shockdat_list[A23fs['id']].P_std)
    axes[ai,aj].errorbar(fsdf['rho'][indfs],fsdf['P'][indfs],fmt='+',color='grey',label='M80 FS excluded')
    indhigh = np.where(fsdf['P'][indfs] > 74)[0]
    axes[ai,aj].errorbar(fsdf['rho'][indfs][indhigh],fsdf['P'][indfs][indhigh],fmt='b+',label='M80 FS included')

    #axes[ai,aj].plot(fitfslrarr,fitfslparr,'b-',linewidth=2,label='cubic fit liquid')
    #axes[ai,aj].plot(fitfssrarr,fitfssparr,'b-',label='linear fit solid')
    ind = np.where(fs_hug_RTK19.P > 74.)[0]
    axes[ai,aj].plot(fs_hug_RTK19.rho[ind],fs_hug_RTK19.P[ind],'b--',label='R19 FS',linewidth=1)
    axes[ai,aj].errorbar(shockdat_list[R19fs['id']].rho,shockdat_list[R19fs['id']].P,fmt='bd',\
             xerr=shockdat_list[R19fs['id']].rho_std,yerr=shockdat_list[R19fs['id']].P_std,label=R19fs['lab'])

axes[ai,aj].set_ylim(0,300)
axes[ai,aj].set_xlim(3.,6.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-5)





#------------------------------
# entropy-temperature
ai=3
aj=0
axes[ai,aj].plot(NewEOS.vc.Sl*1.E3,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
axes[ai,aj].plot(NewEOS.vc.Sv*1.E3,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
axes[ai,aj].plot(NewEOS.cp.S*1.E3,NewEOS.cp.T,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Sl*1.E3,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
    axes[ai,aj].plot(NewEOS.mc.Ss*1.E3,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.S*1.E3,NewEOS.hug.T,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.S[0]*1.E3,NewEOS.hug.T[0],'o',color=stylestphug[1])


axes[ai,aj].plot(MeloshEOS.vc.Sl*1.E3,MeloshEOS.vc.T,styleanphasemelosh[0],color=styleanphasemelosh[1],linewidth=styleanphasemelosh[2],label=styleanphasemelosh[3])
axes[ai,aj].plot(MeloshEOS.vc.Sv*1.E3,MeloshEOS.vc.T,styleanphasemelosh[0],color=styleanphasemelosh[1],linewidth=styleanphasemelosh[2])
axes[ai,aj].plot(MeloshEOS.cp.S*1.E3,MeloshEOS.cp.T,'o',color=stylephasept[1],markersize=stylephasept[2])


if 1:
    #plot Connolly calcs
    b = np.argsort(C16highs)
    axes[ai,aj].plot(C16highs[b]/1.e3,C16hight[b],'-',color='black',label='C16 HSC-vWRK vapor curve')
    axes[ai,aj].plot(5.382,6074,'o',color='black',markersize=5)
    axes[ai,aj].plot(C16lows/1.e3,C16lowt,'--',color='black',label='C16 HSC-CSRK vapor curve')
    axes[ai,aj].plot(5.322,5548,'o',color='black',markersize=5)
    #plt.plot(C16sils,C16silt,'-',color='orangered',label='C16 SiO(l)-vWRK')
    #plt.plot(5403,5675,'o',color='orangered')
    # plot boiling point at 1 bar
    # note Kraus S_vap too small
    axes[ai,aj].plot([3.554,7.460],[3130,3130],'kD',label='C16 1-bar boiling pt.')

axes[ai,aj].set_ylim(0,1.5*NewEOS.cp.T)
axes[ai,aj].set_xlim(0,14)
axes[ai,aj].set_xlabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_ylabel('Temperature [K]',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].legend(fontsize=sfont-7)


#-------
# pressure - temperature melt curve
ai=3
aj=1
#axes[ai,aj].plot(NewEOS.vc.Pl,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
#axes[ai,aj].plot(NewEOS.vc.Pv,NewEOS.vc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])
#axes[ai,aj].plot(NewEOS.cp.P,NewEOS.cp.T,'o',color=stylephasept[1],markersize=stylephasept[2])
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Pl,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2],label=styleanphase[3])
    axes[ai,aj].plot(NewEOS.mc.Ps,NewEOS.mc.T,styleanphase[0],color=styleanphase[1],linewidth=styleanphase[2])

axes[ai,aj].plot(NewEOS.hug.P,NewEOS.hug.T,stylestphug[0],linewidth=stylestphug[2],color=stylestphug[1],label=stylestphug[3])
axes[ai,aj].plot(NewEOS.hug.P[0],NewEOS.hug.T[0],'o',color=stylestphug[1])
axes[ai,aj].plot(NewEOS.hugo.P,NewEOS.hugo.T,'--',linewidth=styleanFS[2],color=styleanFS[1],label=styleanFS[3])

if 0:
    ind = np.where(stsmcsmooth.P <=Pt3)[0]
    axes[ai,aj].plot(stsmcsmooth.P[ind],stsmcsmooth.T[ind],'m-',linewidth=3,label='Proposed melt curve')
    ind = np.where(stsmcsmooth.P >=Pt3)[0]
    axes[ai,aj].plot(stsmcsmooth.P[ind],stsmcsmooth.T[ind],'m:',linewidth=3)#,label='Proposed melt curve')
    # Millot stv melt curve crossing
    axes[ai,aj].fill_between([470.,530.],[8000.,8000.],[8600.,8600.],alpha=0.5,color='blue',label='MM15')

axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].set_xlabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)
axes[ai,aj].set_ylim(0,9000)
axes[ai,aj].set_xlim(-5,600)
axes[ai,aj].legend(fontsize=sfont-5)




plt.savefig(ORIGMODELNAME+'-vs-Melosh-2007.pdf')

Requested ANEOS Hugoniot initial state from (THUG,RHUG) (red curves in plots below):
Initial Pressure (GPa): 0.0001
Initial Temperature (K): 298.0
Initial Sp. Entropy (MJ/K/kg): 0.00017886
Initial Sp. Internal Energy (MJ/kg): 0.043369
Initial Density (g/cm3): 2.2

# check the gridding of the table around the vapor curve

# Plot aneos vapor curves
sfont = 15
fig, axes = plt.subplots(nrows=5, ncols=2, figsize=(14,26))
plt.subplots_adjust(wspace=0.4)
tcolmin=0.
tcolmax=1.2*NewEOS.cp.T
tpltmax=1.2*NewEOS.cp.T

tskip=1
#-------------
ai=0
aj=0
axes[ai,aj].set_xlabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_ylabel('Temperature [K]',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

axes[ai,aj].set_ylim(0,tpltmax)
axes[ai,aj].set_xlim(0,25)

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.S[it,:]*1.e3,np.zeros(NewEOS.ND)+NewEOS.T[it],'o',markersize=2,color=col)

axes[ai,aj].plot(NewEOS.vc.Sl*1.e3,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.vc.Sv*1.e3,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.cp.S*1.e3,NewEOS.cp.T,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Sl*1.e3,NewEOS.mc.T,color="black")
    axes[ai,aj].plot(NewEOS.mc.Ss*1.e3,NewEOS.mc.T,color="black")


#-------------
ai=0
aj=1
for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.P[it,:],'o',markersize=2,color=col)
axes[ai,aj].loglog(NewEOS.vc.rl,NewEOS.vc.Pl,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.Pv,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.P,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(1.e-25,1000)
axes[ai,aj].set_xlim(1.e-20,10.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


#-------------
ai=1
aj=0


for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.S[it,:]*1.E3,NewEOS.P[it,:],'o',markersize=2,color=col)
axes[ai,aj].semilogy(NewEOS.vc.Sl*1.e3,NewEOS.vc.Pl,color="black")
axes[ai,aj].plot(NewEOS.vc.Sv*1.e3,NewEOS.vc.Pv,color="black")
axes[ai,aj].plot(NewEOS.cp.S*1.e3,NewEOS.cp.P,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Sl*1.e3,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.Ss*1.e3,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(1.e-10,1000)
axes[ai,aj].set_xlim(0,15)
axes[ai,aj].set_xlabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-------------
ai=1
aj=1


for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.U[it,:],NewEOS.P[it,:],'o',markersize=2,color=col)

axes[ai,aj].semilogy(NewEOS.vc.Ul,NewEOS.vc.Pl,color="black")
axes[ai,aj].plot(NewEOS.vc.Uv,NewEOS.vc.Pv,color="black")
axes[ai,aj].plot(NewEOS.cp.U,NewEOS.cp.P,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Ul,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.Us,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(1.e-10,1000)
axes[ai,aj].set_xlim(0,25)
axes[ai,aj].set_xlabel('Specific Energy (MJ/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)



#-------------
ai=2
aj=0

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,np.zeros(NewEOS.ND)+NewEOS.T[it],'o',markersize=2,color=col)
axes[ai,aj].plot(NewEOS.vc.rl,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.T,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.T,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.T,color="black")

axes[ai,aj].set_ylim(0.,20000)
axes[ai,aj].set_xlim(0,30.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


#-------------
ai=2
aj=1

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].loglog(NewEOS.rho,np.zeros(NewEOS.ND)+NewEOS.T[it],'o',markersize=2,color=col)

axes[ai,aj].semilogx(NewEOS.vc.rl,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.T,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.T,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.T,color="black")

axes[ai,aj].set_ylim(10.,1.E8)
axes[ai,aj].set_xlim(1.E-10,100.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-------------
ai=3
aj=0

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.S[it,:]*1000.,'o',markersize=2,color=col)

axes[ai,aj].plot(NewEOS.vc.rl,NewEOS.vc.Sl*1.E3,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.Sv*1.E3,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.S*1.E3,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Sl*1.E3,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ss*1.E3,color="black")


axes[ai,aj].set_ylim(0.,10)
axes[ai,aj].set_xlim(0.,15.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-------------
ai=3
aj=1

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].semilogx(NewEOS.rho,NewEOS.S[it,:]*1000.,'o',markersize=2,color=col)

axes[ai,aj].semilogx(NewEOS.vc.rl,NewEOS.vc.Sl*1.E3,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.Sv*1.E3,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.S*1.E3,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Sl*1.E3,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ss*1.E3,color="black")

axes[ai,aj].set_ylim(0.,20)
axes[ai,aj].set_xlim(1.E-15,30.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)




#-------------
ai=4
aj=1

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].semilogx(NewEOS.rho,NewEOS.S[it,:]*1000.,'o',markersize=2,color=col)

axes[ai,aj].semilogx(NewEOS.vc.rl,NewEOS.vc.Sl*1.E3,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.Sv*1.E3,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.S*1.E3,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Sl*1.E3,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ss*1.E3,color="black")

axes[ai,aj].set_ylim(0.,30)
axes[ai,aj].set_xlim(1.E-25,30.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-------------
ai=4
aj=0

tcolmax=40000.
for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].semilogx(NewEOS.rho,NewEOS.S[it,:]*1000.,'o',markersize=2,color=col)

axes[ai,aj].semilogx(NewEOS.vc.rl,NewEOS.vc.Sl*1.E3,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.Sv*1.E3,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.S*1.E3,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Sl*1.E3,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ss*1.E3,color="black")

axes[ai,aj].set_ylim(0.,200)
axes[ai,aj].set_xlim(1.E-25,30.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#axes[0,0].legend(fontsize=sfont-5)
#axes[0,1].legend(fontsize=sfont-5)
#axes[1,0].legend(fontsize=sfont-5)
#axes[1,1].legend(fontsize=sfont-5)

# check the ANEOS quality flag for errors in the EOS surface
uniqkarr = np.unique(NewEOS.KPA)
print("Unique values in the KPA flag array:", uniqkarr)
print("Value: # elements with this value")
#get_indexes = lambda x, xs: [i for (y, i) in zip(xs, range(len(xs))) if x == y]

for i in uniqkarr:
    ix = get_indexes(i,np.reshape(NewEOS.KPA,NewEOS.ND*NewEOS.NT))
    print(i,': ',len(ix))


    
# Plot aneos vapor curves
sfont = 15
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(14,14))
plt.subplots_adjust(wspace=0.4)
kcolmin=1
kcolmax=7

tskip=10
#-------------


#-------------
ai=0
aj=0
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


for itt in range(1,NewEOS.NT,tskip):
    for i in uniqkarr:
        ix = get_indexes(i,NewEOS.KPA[itt,:])
        icolscale = (i-kcolmin)/(kcolmax-kcolmin)
        col=local_cmaps.parula(icolscale)
        if i == 5:
            axes[ai,aj].plot(NewEOS.rho[ix],np.full(len(ix),NewEOS.T[itt]),'o',markersize=5,color='red')
#            print(i,andarr[ix],antarr[itt])
        else:
            axes[ai,aj].plot(NewEOS.rho[ix],np.full(len(ix),NewEOS.T[itt]),'o',markersize=2,color=col)

axes[ai,aj].plot(NewEOS.vc.rl,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.T,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.T,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.T,color="black")


axes[ai,aj].set_ylim(0.,20000)
axes[ai,aj].set_xlim(0,30.)


#-------------
ai=0
aj=1

for itt in range(1,NewEOS.NT,tskip):
    for i in uniqkarr:
        ix = get_indexes(i,NewEOS.KPA[itt,:])
        icolscale = (i-kcolmin)/(kcolmax-kcolmin)
        col=local_cmaps.parula(icolscale)
        if i == 5:
            axes[ai,aj].loglog(NewEOS.rho[ix],np.full(len(ix),NewEOS.T[itt]),'o',markersize=5,color='red')
#            print(i,andarr[ix],antarr[itt])
        else:
            axes[ai,aj].loglog(NewEOS.rho[ix],np.full(len(ix),NewEOS.T[itt]),'o',markersize=2,color=col)

axes[ai,aj].loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.T,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.T,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.T,color="black")

axes[ai,aj].set_ylim(10.,1.E8)
axes[ai,aj].set_xlim(1.E-22,100.)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


#-------------
ai=1
aj=0
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


for itt in range(1,NewEOS.NT,tskip):
    for i in uniqkarr:
        ix = get_indexes(i,NewEOS.KPA[itt,:])
        icolscale = (i-kcolmin)/(kcolmax-kcolmin)
        col=local_cmaps.parula(icolscale)
        if i == 5:
            axes[ai,aj].semilogy(NewEOS.rho[ix],NewEOS.P[itt,ix],'o',markersize=5,color='red')
#            print(i,andarr[ix],antarr[itt])
        else:
            axes[ai,aj].semilogy(NewEOS.rho[ix],NewEOS.P[itt,ix],'o',markersize=2,color=col)

axes[ai,aj].semilogy(NewEOS.vc.rl,NewEOS.vc.Pl,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.Pv,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.P,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(1.e-6,2000)
axes[ai,aj].set_xlim(0,13.)


#-------------
ai=1
aj=1
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Temperature (K)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


for itt in range(1,NewEOS.NT,tskip):
    for i in uniqkarr:
        ix = get_indexes(i,NewEOS.KPA[itt,:])
        icolscale = (i-kcolmin)/(kcolmax-kcolmin)
        col=local_cmaps.parula(icolscale)
        if i == 5:
            axes[ai,aj].plot(NewEOS.rho[ix],np.full(len(ix),NewEOS.T[itt]),'o',markersize=5,color='red')
#            print(i,andarr[ix],antarr[itt])
        else:
            axes[ai,aj].plot(NewEOS.rho[ix],np.full(len(ix),NewEOS.T[itt]),'o',markersize=2,color=col)
    tmp = np.where(NewEOS.P[itt,:] < 0.)
    if (len(tmp[0]) > 0):
        axes[ai,aj].plot(NewEOS.rho[tmp[0]],np.full(len(tmp[0]),NewEOS.T[itt]),'o',markersize=2,color='purple')
    tmp = np.where(NewEOS.P[itt,:]*1.e30 <= 2.)
    if (len(tmp[0]) > 0):
        axes[ai,aj].plot(NewEOS.rho[tmp[0]],np.full(len(tmp[0]),NewEOS.T[itt]),'o',markersize=2,color='hotpink')

axes[ai,aj].loglog(NewEOS.vc.rl,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.vc.rv,NewEOS.vc.T,color="black")
axes[ai,aj].plot(NewEOS.cp.rho,NewEOS.cp.T,'o',color="black")
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.T,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.T,color="black")
    
axes[ai,aj].set_ylim(10.,1.E8)
axes[ai,aj].set_xlim(1.E-22,100.)
#fig.show()

Unique values in the KPA flag array: [2. 4. 5. 6.]
Value: # elements with this value
2.0 :  104934
4.0 :  213143
5.0 :  935
6.0 :  246665

(1e-22, 100.0)

# CHECK MELT CURVE GRIDDING
sfont = 15
fig, axes = plt.subplots(nrows=5, ncols=2, figsize=(14,26))
plt.subplots_adjust(wspace=0.5)
tcolmin=0.
tcolmax=1.5*NewEOS.cp.T
tpltmax=1.5*NewEOS.cp.T

tskip=10
#-------------
ai=0
aj=0
axes[ai,aj].plot([],[],' ',label='Low pressure melt curve')

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(.1,60)
axes[ai,aj].set_xlim(2.4,4.5)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-------------
ai=0
aj=1
axes[ai,aj].plot([],[],' ',label='High pressure melt curve')

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(40.,150)
axes[ai,aj].set_xlim(3.9,5.4)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


#-------------
ai=1
aj=0
axes[ai,aj].plot([],[],' ',label='Higher pressure melt curve')

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,color="black")


axes[ai,aj].set_ylim(180.,310)
axes[ai,aj].set_xlim(5.5,6.4)
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


#-------------
ai=1
aj=1
axes[ai,aj].plot([],[],' ',label='Higher pressure melt curve')

axes[ai,aj].set_xlabel('Specific Energy (MJ/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.U[it,:],NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Ul,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.Us,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(180.,310)
axes[ai,aj].set_xlim(7,20)

#-------------

ai=2
aj=0
axes[ai,aj].plot([],[],' ',label='Low pressure melt curve')

axes[ai,aj].set_xlabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.S[it,:]*1.e3,NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Sl*1.e3,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.Ss*1.e3,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(.1,50)
axes[ai,aj].set_xlim(2.6,4.)


#-------------

ai=2
aj=1
axes[ai,aj].plot([],[],' ',label='High pressure melt curve')

axes[ai,aj].set_xlabel('Specific Entropy (kJ/K/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.S[it,:]*1.e3,NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Sl*1.e3,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.Ss*1.e3,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(50.,150)
axes[ai,aj].set_xlim(2.6,4.2)


#-------------
ai=3
aj=0
axes[ai,aj].plot([],[],' ',label='Low pressure melt curve')

axes[ai,aj].set_xlabel('Specific Energy (MJ/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.U[it,:],NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Ul,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.Us,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_xlim(1,7)
axes[ai,aj].set_ylim(.1,50)

#-------------
ai=3
aj=1
axes[ai,aj].plot([],[],' ',label='High pressure melt curve')

axes[ai,aj].set_xlabel('Specific Energy (MJ/kg)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)


for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.U[it,:],NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.Ul,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.Us,NewEOS.mc.Ps,color="black")
axes[ai,aj].set_ylim(50.,150)
axes[ai,aj].set_xlim(1,15)



#-------------
ai=4
aj=0
for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.P[it,:],'o',markersize=2,color=col,linestyle='-')

if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,color="black")

axes[ai,aj].set_ylim(300,700)
axes[ai,aj].set_xlim(6.5,8.)

axes[ai,aj].plot([],[],' ',label='Very high pressure melt curve.')
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#-------------
ai=4
aj=1
for it in range(1,NewEOS.NT,tskip):
    icolscale = (NewEOS.T[it]-tcolmin)/(tcolmax-tcolmin)
    col=local_cmaps.parula(icolscale)
    axes[ai,aj].plot(NewEOS.rho,NewEOS.P[it,:],'o',markersize=2,linestyle='-',color=col)
if NewEOS.tp.T > 0:
    axes[ai,aj].plot(NewEOS.mc.rl,NewEOS.mc.Pl,color="black")
    axes[ai,aj].plot(NewEOS.mc.rs,NewEOS.mc.Ps,color="black")
axes[ai,aj].set_ylim(700,1200)
axes[ai,aj].set_xlim(8.25,9.5)

axes[ai,aj].plot([],[],' ',label='Very high pressures')
axes[ai,aj].set_xlabel('Density (g/cm$^3$)',size=sfont)
axes[ai,aj].set_ylabel('Pressure (GPa)',size=sfont)
axes[ai,aj].tick_params(labelsize=sfont)

#------
axes[0,0].legend(fontsize=sfont-5)
axes[0,1].legend(fontsize=sfont-5)
axes[1,0].legend(fontsize=sfont-5)
axes[1,1].legend(fontsize=sfont-5)
axes[2,0].legend(fontsize=sfont-5)
axes[2,1].legend(fontsize=sfont-5)
axes[3,0].legend(fontsize=sfont-5)
axes[3,1].legend(fontsize=sfont-5)
axes[4,0].legend(fontsize=sfont-5)
axes[4,1].legend(fontsize=sfont-5)

<matplotlib.legend.Legend at 0x336595e90>

ANEOS SILICA MODEL DOCUMENTATION NOTEBOOK¶

ANEOS FUSED SILICA v0.2¶

USER INPUTS FOR SESAME AND GADGET2 TABLE CONSTRUCTION¶

Read in silica data¶

Color mesh plots of the SESAME Rho-T table¶

Model Data Quality (MDQ) Flag¶

ANEOS Model Compared to Experimental and Ab Initio Data¶

Model Critical Point and Shock-Induced Phase Changes¶

Model Phase Boundaries¶

ANEOS model comparisons to experimental data¶

ANEOS Gruneisen parameter and Theta parameter for thermal model¶

Specific Heat Capacity and Sounds Speeds¶

Compare Silica EOS Models¶

Check the grid spacing compared to the phase boundaries¶

ANEOS KPA FLAG¶

Check the Melt Curve in Tabulated EOS¶

References¶

ANEOS references¶

ANEOS Input Parameters¶

Development Notes¶

¶